Changeset 4897
- Timestamp:
- Dec 27, 2007, 6:35:08 PM (17 years ago)
- Location:
- anuga_core/source/anuga
- Files:
-
- 5 edited
Legend:
- Unmodified
- Added
- Removed
-
anuga_core/source/anuga/abstract_2d_finite_volumes/domain.py
r4836 r4897 380 380 Depth = domain.create_quantity_from_expression('stage-elevation') 381 381 382 exp = '(xmomentum*xmomentum + ymomentum*ymomentum)**0.5' )382 exp = '(xmomentum*xmomentum + ymomentum*ymomentum)**0.5' 383 383 Absolute_momentum = domain.create_quantity_from_expression(exp) 384 384 -
anuga_core/source/anuga/abstract_2d_finite_volumes/quantity.py
r4886 r4897 250 250 polygon = None, 251 251 indices = None, 252 smooth = False, 252 253 verbose = False, 253 254 use_cache = False): … … 389 390 self.set_values_from_constant(numeric, 390 391 location, indices, verbose) 391 392 393 392 394 self.extrapolate_first_order() 395 396 if smooth: 397 self.smooth_vertex_values() 398 399 393 400 return 394 401 … … 1181 1188 1182 1189 1183 def smooth_vertex_values(self, value_array='field_values', 1184 precision = None): 1185 """ Smooths field_values or conserved_quantities data. 1186 TODO: be able to smooth individual fields 1187 NOTE: This function does not have a test. 1188 FIXME: NOT DONE - do we need it? 1189 FIXME: this function isn't called by anything. 1190 Maybe it should be removed..-DSG 1191 """ 1192 1193 from Numeric import concatenate, zeros, Float, Int, array, reshape 1194 1195 1196 A,V = self.get_vertex_values(xy=False, 1197 value_array=value_array, 1198 smooth = True, 1199 precision = precision) 1200 1201 #Set some field values 1202 for volume in self: 1203 for i,v in enumerate(volume.vertices): 1204 if value_array == 'field_values': 1205 volume.set_field_values('vertex', i, A[v,:]) 1206 elif value_array == 'conserved_quantities': 1207 volume.set_conserved_quantities('vertex', i, A[v,:]) 1208 1209 if value_array == 'field_values': 1210 self.precompute() 1211 elif value_array == 'conserved_quantities': 1212 Volume.interpolate_conserved_quantities() 1190 def smooth_vertex_values(self): 1191 """ Smooths vertex values. 1192 """ 1193 1194 A,V = self.get_vertex_values(xy=False, smooth = True) 1195 self.set_vertex_values(A) 1213 1196 1214 1197 … … 1254 1237 if smooth is None: 1255 1238 # Take default from domain 1256 smooth = self.domain.smooth 1239 try: 1240 smooth = self.domain.smooth 1241 except: 1242 smooth = False 1257 1243 1258 1244 if precision is None: … … 1392 1378 return update(self, timestep) 1393 1379 1394 1395 1380 def compute_gradients(self): 1396 1381 #Call correct module function … … 1398 1383 return compute_gradients(self) 1399 1384 1400 1401 1385 def limit(self): 1402 1386 #Call correct module depending on whether … … 1404 1388 limit_old(self) 1405 1389 1406 1407 def limit_by_vertex(self): 1390 def limit_vertices_by_all_neighbours(self): 1408 1391 #Call correct module function 1409 1392 #(either from this module or C-extension) 1410 limit_by_vertex(self) 1411 1412 1413 def limit_by_edge(self): 1393 limit_vertices_by_all_neighbours(self) 1394 1395 def limit_edges_by_all_neighbours(self): 1414 1396 #Call correct module function 1415 1397 #(either from this module or C-extension) 1416 limit_by_edge(self) 1398 limit_edges_by_all_neighbours(self) 1399 1400 def limit_edges_by_neighbour(self): 1401 #Call correct module function 1402 #(either from this module or C-extension) 1403 limit_edges_by_neighbour(self) 1417 1404 1418 1405 def extrapolate_second_order(self): … … 1420 1407 #(either from this module or C-extension) 1421 1408 extrapolate_second_order(self) 1422 1423 1409 1424 1410 def backup_centroid_values(self): … … 1444 1430 compute_gradients,\ 1445 1431 limit_old,\ 1446 limit_by_vertex,\ 1447 limit_by_edge,\ 1432 limit_vertices_by_all_neighbours,\ 1433 limit_edges_by_all_neighbours,\ 1434 limit_edges_by_neighbour,\ 1448 1435 extrapolate_second_order,\ 1449 1436 interpolate_from_vertices_to_edges,\ -
anuga_core/source/anuga/abstract_2d_finite_volumes/quantity_ext.c
r4886 r4897 137 137 138 138 139 int _limit_ by_vertex(int N, double beta,139 int _limit_vertices_by_all_neighbours(int N, double beta, 140 140 double* centroid_values, 141 141 double* vertex_values, … … 195 195 } 196 196 197 int _limit_ by_edge(int N, double beta,197 int _limit_edges_by_all_neighbours(int N, double beta, 198 198 double* centroid_values, 199 199 double* vertex_values, … … 240 240 241 241 //Update vertex and edge values using phi limiter 242 edge_values[k3+0] = qc + phi*dqa[0]; 243 edge_values[k3+1] = qc + phi*dqa[1]; 244 edge_values[k3+2] = qc + phi*dqa[2]; 245 246 vertex_values[k3+0] = edge_values[k3+1] + edge_values[k3+2] - edge_values[k3+0]; 247 vertex_values[k3+1] = edge_values[k3+2] + edge_values[k3+0] - edge_values[k3+1]; 248 vertex_values[k3+2] = edge_values[k3+0] + edge_values[k3+1] - edge_values[k3+2]; 249 250 } 251 252 return 0; 253 } 254 255 256 int _limit_edges_by_neighbour(int N, double beta, 257 double* centroid_values, 258 double* vertex_values, 259 double* edge_values, 260 long* neighbours) { 261 262 int i, k, k2, k3, k6; 263 long n; 264 double qmin, qmax, qn, qc; 265 double dq, dqa[3], phi, r; 266 267 for (k=0; k<N; k++){ 268 k6 = 6*k; 269 k3 = 3*k; 270 k2 = 2*k; 271 272 qc = centroid_values[k]; 273 phi = 1.0; 274 275 for (i=0; i<3; i++) { 276 dq = edge_values[k3+i] - qc; //Delta between edge and centroid values 277 dqa[i] = dq; //Save dq for use in updating vertex values 278 279 n = neighbours[k3+i]; 280 if (n >= 0) { 281 qn = centroid_values[n]; //Neighbour's centroid value 282 283 qmin = min(qc, qn); 284 qmax = max(qc, qn); 285 286 r = 1.0; 287 288 if (dq > 0.0) r = (qmax - qc)/dq; 289 if (dq < 0.0) r = (qmin - qc)/dq; 290 291 phi = min( min(r*beta, 1.0), phi); 292 } 293 } 294 295 296 //Update edge and vertex values using phi limiter 242 297 edge_values[k3+0] = qc + phi*dqa[0]; 243 298 edge_values[k3+1] = qc + phi*dqa[1]; … … 974 1029 975 1030 976 PyObject *limit_ by_vertex(PyObject *self, PyObject *args) {1031 PyObject *limit_vertices_by_all_neighbours(PyObject *self, PyObject *args) { 977 1032 //Limit slopes for each volume to eliminate artificial variance 978 1033 //introduced by e.g. second order extrapolator … … 1033 1088 N = centroid_values -> dimensions[0]; 1034 1089 1035 err = _limit_ by_vertex(N, beta_w,1036 1037 1038 1039 1040 1090 err = _limit_vertices_by_all_neighbours(N, beta_w, 1091 (double*) centroid_values -> data, 1092 (double*) vertex_values -> data, 1093 (double*) edge_values -> data, 1094 (long*) neighbours -> data); 1095 1041 1096 if (err != 0) { 1042 1097 PyErr_SetString(PyExc_RuntimeError, … … 1059 1114 1060 1115 1061 PyObject *limit_ by_edge(PyObject *self, PyObject *args) {1116 PyObject *limit_edges_by_all_neighbours(PyObject *self, PyObject *args) { 1062 1117 //Limit slopes for each volume to eliminate artificial variance 1063 1118 //introduced by e.g. second order extrapolator … … 1086 1141 if (!PyArg_ParseTuple(args, "O", &quantity)) { 1087 1142 PyErr_SetString(PyExc_RuntimeError, 1088 "quantity_ext.c: limit_ by_edgecould not parse input");1143 "quantity_ext.c: limit_edges_by_all_neighbours could not parse input"); 1089 1144 return NULL; 1090 1145 } … … 1093 1148 if (!domain) { 1094 1149 PyErr_SetString(PyExc_RuntimeError, 1095 "quantity_ext.c: limit_ by_edgecould not obtain domain object from quantity");1150 "quantity_ext.c: limit_edges_by_all_neighbours could not obtain domain object from quantity"); 1096 1151 1097 1152 return NULL; … … 1102 1157 if (!Tmp) { 1103 1158 PyErr_SetString(PyExc_RuntimeError, 1104 "quantity_ext.c: limit_ by_edgecould not obtain beta_w object from domain");1159 "quantity_ext.c: limit_edges_by_all_neighbours could not obtain beta_w object from domain"); 1105 1160 1106 1161 return NULL; … … 1118 1173 N = centroid_values -> dimensions[0]; 1119 1174 1120 err = _limit_ by_edge(N, beta_w,1121 1122 1123 1124 1175 err = _limit_edges_by_all_neighbours(N, beta_w, 1176 (double*) centroid_values -> data, 1177 (double*) vertex_values -> data, 1178 (double*) edge_values -> data, 1179 (long*) neighbours -> data); 1125 1180 1126 1181 if (err != 0) { … … 1143 1198 1144 1199 1200 PyObject *limit_edges_by_neighbour(PyObject *self, PyObject *args) { 1201 //Limit slopes for each volume to eliminate artificial variance 1202 //introduced by e.g. second order extrapolator 1203 1204 //This is an unsophisticated limiter as it does not take into 1205 //account dependencies among quantities. 1206 1207 //precondition: 1208 // vertex values are estimated from gradient 1209 //postcondition: 1210 // vertex and edge values are updated 1211 // 1212 1213 PyObject *quantity, *domain, *Tmp; 1214 PyArrayObject 1215 *vertex_values, //Conserved quantities at vertices 1216 *centroid_values, //Conserved quantities at centroids 1217 *edge_values, //Conserved quantities at edges 1218 *neighbours; 1219 1220 double beta_w; //Safety factor 1221 int N, err; 1222 1223 1224 // Convert Python arguments to C 1225 if (!PyArg_ParseTuple(args, "O", &quantity)) { 1226 PyErr_SetString(PyExc_RuntimeError, 1227 "quantity_ext.c: limit_edges_by_neighbour could not parse input"); 1228 return NULL; 1229 } 1230 1231 domain = PyObject_GetAttrString(quantity, "domain"); 1232 if (!domain) { 1233 PyErr_SetString(PyExc_RuntimeError, 1234 "quantity_ext.c: limit_edges_by_neighbour could not obtain domain object from quantity"); 1235 1236 return NULL; 1237 } 1238 1239 // Get safety factor beta_w 1240 Tmp = PyObject_GetAttrString(domain, "beta_w"); 1241 if (!Tmp) { 1242 PyErr_SetString(PyExc_RuntimeError, 1243 "quantity_ext.c: limit_by_vertex could not obtain beta_w object from domain"); 1244 1245 return NULL; 1246 } 1247 1248 1249 // Get pertinent variables 1250 neighbours = get_consecutive_array(domain, "neighbours"); 1251 centroid_values = get_consecutive_array(quantity, "centroid_values"); 1252 vertex_values = get_consecutive_array(quantity, "vertex_values"); 1253 edge_values = get_consecutive_array(quantity, "edge_values"); 1254 beta_w = PyFloat_AsDouble(Tmp); 1255 1256 1257 N = centroid_values -> dimensions[0]; 1258 1259 err = _limit_edges_by_neighbour(N, beta_w, 1260 (double*) centroid_values -> data, 1261 (double*) vertex_values -> data, 1262 (double*) edge_values -> data, 1263 (long*) neighbours -> data); 1264 1265 if (err != 0) { 1266 PyErr_SetString(PyExc_RuntimeError, 1267 "Internal function _limit_by_vertex failed"); 1268 return NULL; 1269 } 1270 1271 1272 // Release 1273 Py_DECREF(neighbours); 1274 Py_DECREF(centroid_values); 1275 Py_DECREF(vertex_values); 1276 Py_DECREF(edge_values); 1277 Py_DECREF(Tmp); 1278 1279 1280 return Py_BuildValue(""); 1281 } 1282 1145 1283 1146 1284 // Method table for python module 1147 1285 static struct PyMethodDef MethodTable[] = { 1148 1286 {"limit_old", limit_old, METH_VARARGS, "Print out"}, 1149 {"limit_by_vertex", limit_by_vertex, METH_VARARGS, "Print out"}, 1150 {"limit_by_edge", limit_by_edge, METH_VARARGS, "Print out"}, 1287 {"limit_vertices_by_all_neighbours", limit_vertices_by_all_neighbours, METH_VARARGS, "Print out"}, 1288 {"limit_edges_by_all_neighbours", limit_edges_by_all_neighbours, METH_VARARGS, "Print out"}, 1289 {"limit_edges_by_neighbour", limit_edges_by_neighbour, METH_VARARGS, "Print out"}, 1151 1290 {"update", update, METH_VARARGS, "Print out"}, 1152 1291 {"backup_centroid_values", backup_centroid_values, METH_VARARGS, "Print out"}, -
anuga_core/source/anuga/abstract_2d_finite_volumes/test_quantity.py
r4886 r4897 1445 1445 1446 1446 1447 def test_ vertex_limiter(self):1447 def test_limit_vertices_by_all_neighbours(self): 1448 1448 quantity = Conserved_quantity(self.mesh4) 1449 1449 … … 1453 1453 1454 1454 #Limit 1455 quantity.limit_ by_vertex()1455 quantity.limit_vertices_by_all_neighbours() 1456 1456 1457 1457 #Assert that central triangle is limited by neighbours … … 1475 1475 1476 1476 1477 def test_ edge_limiter(self):1477 def test_limit_edges_by_all_neighbours(self): 1478 1478 quantity = Conserved_quantity(self.mesh4) 1479 1479 … … 1483 1483 1484 1484 #Limit 1485 quantity.limit_ by_edge()1485 quantity.limit_edges_by_all_neighbours() 1486 1486 1487 1487 #Assert that central triangle is limited by neighbours … … 1493 1493 1494 1494 assert quantity.edge_values[1,2] <= quantity.centroid_values[2] 1495 assert quantity.edge_values[1,2] <= quantity.centroid_values[0]1495 assert quantity.edge_values[1,2] >= quantity.centroid_values[0] 1496 1496 1497 1497 … … 1503 1503 sum(quantity.vertex_values[k,:])/3) 1504 1504 1505 1506 def test_limit_edges_by_neighbour(self): 1507 quantity = Conserved_quantity(self.mesh4) 1508 1509 #Create a deliberate overshoot (e.g. from gradient computation) 1510 quantity.set_values([[3,0,3], [2,2,6], [5,3,8], [8,3,5]]) 1511 1512 1513 #Limit 1514 quantity.limit_edges_by_neighbour() 1515 1516 #Assert that central triangle is limited by neighbours 1517 assert quantity.edge_values[1,0] <= quantity.centroid_values[3] 1518 assert quantity.edge_values[1,0] >= quantity.centroid_values[1] 1519 1520 assert quantity.edge_values[1,1] <= quantity.centroid_values[2] 1521 assert quantity.edge_values[1,1] >= quantity.centroid_values[1] 1522 1523 assert quantity.edge_values[1,2] <= quantity.centroid_values[1] 1524 assert quantity.edge_values[1,2] >= quantity.centroid_values[0] 1525 1526 1527 1528 #Assert that quantities are conserved 1529 from Numeric import sum 1530 for k in range(quantity.centroid_values.shape[0]): 1531 assert allclose (quantity.centroid_values[k], 1532 sum(quantity.vertex_values[k,:])/3) 1505 1533 1506 1534 def test_limiter2(self): … … 2207 2235 assert allclose(subset, answer) 2208 2236 2237 def test_smooth_vertex_values(self): 2238 """ 2239 get values based on triangle lists. 2240 """ 2241 from mesh_factory import rectangular 2242 from shallow_water import Domain 2243 from Numeric import zeros, Float 2244 2245 #Create basic mesh 2246 points, vertices, boundary = rectangular(2, 2) 2247 2248 #print "points",points 2249 #print "vertices",vertices 2250 #print "boundary",boundary 2251 2252 #Create shallow water domain 2253 domain = Domain(points, vertices, boundary) 2254 #print "domain.number_of_elements ",domain.number_of_elements 2255 quantity = Quantity(domain,[[0,0,0],[1,1,1],[2,2,2],[3,3,3], 2256 [4,4,4],[5,5,5],[6,6,6],[7,7,7]]) 2257 2258 #print "quantity.get_values(location = 'unique vertices')", \ 2259 # quantity.get_values(location = 'unique vertices') 2260 2261 #print "quantity.get_values(location = 'unique vertices')", \ 2262 # quantity.get_values(indices=[0,1,2,3,4,5,6,7], \ 2263 # location = 'unique vertices') 2264 2265 #print quantity.get_values(location = 'unique vertices') 2266 #print quantity.domain.number_of_triangles_per_node 2267 #print quantity.vertex_values 2268 2269 #answer = [0.5, 2, 3, 3, 3.5, 4, 4, 5, 6.5] 2270 #assert allclose(answer, 2271 # quantity.get_values(location = 'unique vertices')) 2272 2273 quantity.smooth_vertex_values() 2274 2275 #print quantity.vertex_values 2276 2277 2278 answer_vertex_values = [[3,3.5,0.5],[2,0.5,3.5],[3.5,4,2],[3,2,4], 2279 [4,5,3],[3.5,3,5],[5,6.5,3.5],[4,3.5,6.5]] 2280 2281 assert allclose(answer_vertex_values, 2282 quantity.vertex_values) 2283 #print "quantity.centroid_values",quantity.centroid_values 2284 #print "quantity.get_values(location = 'centroids') ",\ 2285 # quantity.get_values(location = 'centroids') 2209 2286 2210 2287 -
anuga_core/source/anuga/visualiser/realtime.py
r3966 r4897 1 from Numeric import Float, zeros 1 from Numeric import Float, zeros, shape 2 2 from Tkinter import Button, E, Tk, W 3 3 from threading import Event … … 45 45 verticies = self.source.get_vertex_coordinates() 46 46 N_vert = len(verticies) 47 47 48 # Also build vert_index - a list of the x & y values of each vertex 48 49 self.vert_index = zeros((N_vert,2), Float) … … 50 51 self.vtk_cells.InsertNextCell(3) 51 52 for v in range(3): 52 self.vert_index[ triangles[n][v]] = verticies[n * 3 + v]53 self.vtk_cells.InsertCellPoint( triangles[n][v])53 self.vert_index[n * 3 + v] = verticies[n * 3 + v] 54 self.vtk_cells.InsertCellPoint(n * 3 + v) 54 55 55 56 def update_height_quantity(self, quantityName, dynamic=True): … … 59 60 vertex_values, _ = self.source.get_quantity(quantityName).get_vertex_values(xy=False, smooth=False) 60 61 61 for n in range(len(triangles)): 62 for v in range(3): 63 #qty_index[triangles[n][v]] = self.source.get_quantity(quantityName).vertex_values[n][v] 64 qty_index[triangles[n][v]] = vertex_values[n * 3 + v] 62 for n in range(N_vert): 63 qty_index[n] = vertex_values[n] 65 64 66 65 points = vtkPoints() … … 124 123 Visualiser.redraw(self) 125 124 126 def update(self ):125 def update(self,pause=False): 127 126 """Sync the visualiser to the domain. Call this in the evolve loop.""" 127 128 128 if self.running: 129 129 self.sync_redrawReady.set() … … 132 132 self.sync_unpaused.wait() 133 133 134 if pause and self.running: 135 if self.sync_unpaused.isSet(): 136 self.sync_unpaused.clear() 137 self.tk_pauseResume.config(text="Resume") 138 139 self.sync_redrawReady.set() 140 self.sync_idle.wait() 141 self.sync_idle.clear() 142 self.sync_unpaused.wait() 143 144 145 134 146 def evolveFinished(self): 135 147 """Stop the visualiser from waiting on signals from the evolve loop.
Note: See TracChangeset
for help on using the changeset viewer.