Changeset 3678
- Timestamp:
- Oct 1, 2006, 6:55:46 PM (18 years ago)
- Files:
-
- 17 edited
Legend:
- Unmodified
- Added
- Removed
-
anuga_core/source/anuga/abstract_2d_finite_volumes/domain.py
r3631 r3678 166 166 self.set_default_order(1) 167 167 #self.default_order = 1 168 #self. order= self.default_order168 #self._order_ = self.default_order 169 169 170 170 self.smallsteps = 0 … … 222 222 223 223 self.default_order = n 224 self. order= self.default_order224 self._order_ = self.default_order 225 225 226 226 … … 695 695 yieldstep = float(yieldstep) 696 696 697 self. order= self.default_order697 self._order_ = self.default_order 698 698 699 699 … … 759 759 self.yieldtime += self.timestep 760 760 self.number_of_steps += 1 761 if self. order== 1:761 if self._order_ == 1: 762 762 self.number_of_first_order_steps += 1 763 763 … … 846 846 self.smallsteps = 0 #Reset 847 847 848 if self. order== 1:848 if self._order_ == 1: 849 849 msg = 'WARNING: Too small timestep %.16f reached '\ 850 850 %timestep … … 857 857 else: 858 858 #Try to overcome situation by switching to 1 order 859 self. order= 1859 self._order_ = 1 860 860 861 861 else: 862 862 self.smallsteps = 0 863 if self. order== 1 and self.default_order == 2:864 self. order= 2863 if self._order_ == 1 and self.default_order == 2: 864 self._order_ = 2 865 865 866 866 … … 933 933 for name in self.conserved_quantities: 934 934 Q = self.quantities[name] 935 if self. order== 1:935 if self._order_ == 1: 936 936 Q.extrapolate_first_order() 937 elif self. order== 2:937 elif self._order_ == 2: 938 938 Q.extrapolate_second_order() 939 939 Q.limit() -
anuga_core/source/anuga/abstract_2d_finite_volumes/mesh_factory.py
r3560 r3678 334 334 """ 335 335 336 from anuga. abstract_2d_finite_volumes.utilimport anglediff336 from anuga.utilities.numerical_tools import anglediff 337 337 from math import pi 338 338 import os.path … … 449 449 450 450 from math import pi 451 from anuga. abstract_2d_finite_volumes.utilimport anglediff451 from anuga.utilities.numerical_tools import anglediff 452 452 453 453 -
anuga_core/source/anuga/abstract_2d_finite_volumes/quantity.py
r3560 r3678 1234 1234 #Get index of the one neighbour 1235 1235 for k0 in surrogate_neighbours[k,:]: 1236 1236 if k0 != k: break 1237 1237 assert k0 != k 1238 1238 … … 1285 1285 1286 1286 for k in range(N): 1287 qmax[k] = qmin[k] = qc[k] 1287 qmax[k] = qc[k] 1288 qmin[k] = qc[k] 1288 1289 for i in range(3): 1289 1290 n = quantity.domain.neighbours[k,i] … … 1293 1294 qmin[k] = min(qmin[k], qn) 1294 1295 qmax[k] = max(qmax[k], qn) 1296 qmax[k] = min(qmax[k], 2.0*qc[k]) 1297 qmin[k] = max(qmin[k], 0.5*qc[k]) 1295 1298 1296 1299 … … 1326 1329 #Replace python version with c implementations 1327 1330 1328 from quantity_ext import limit, compute_gradients,\1331 from quantity_ext import compute_gradients, limit,\ 1329 1332 extrapolate_second_order, interpolate_from_vertices_to_edges, update -
anuga_core/source/anuga/abstract_2d_finite_volumes/quantity_ext.c
r2676 r3678 26 26 long* surrogate_neighbours, 27 27 double* a, 28 double* b) 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60 61 62 63 64 65 //#Get index of the one neighbour66 67 68 k0 = surrogate_neighbours[index3 + i];69 i++;70 71 72 73 74 75 76 77 78 79 80 81 82 83 84 85 86 87 88 89 90 91 92 93 94 95 96 97 98 28 double* b){ 29 30 int i, k, k0, k1, k2, index3; 31 double x0, x1, x2, y0, y1, y2, q0, q1, q2; //, det; 32 33 34 for (k=0; k<N; k++) { 35 index3 = 3*k; 36 37 if (number_of_boundaries[k] < 2) { 38 //Two or three true neighbours 39 40 //Get indices of neighbours (or self when used as surrogate) 41 //k0, k1, k2 = surrogate_neighbours[k,:] 42 43 k0 = surrogate_neighbours[index3 + 0]; 44 k1 = surrogate_neighbours[index3 + 1]; 45 k2 = surrogate_neighbours[index3 + 2]; 46 47 48 if (k0 == k1 || k1 == k2) return -1; 49 50 //Get data 51 q0 = centroid_values[k0]; 52 q1 = centroid_values[k1]; 53 q2 = centroid_values[k2]; 54 55 x0 = centroids[k0*2]; y0 = centroids[k0*2+1]; 56 x1 = centroids[k1*2]; y1 = centroids[k1*2+1]; 57 x2 = centroids[k2*2]; y2 = centroids[k2*2+1]; 58 59 //Gradient 60 _gradient(x0, y0, x1, y1, x2, y2, q0, q1, q2, &a[k], &b[k]); 61 62 } else if (number_of_boundaries[k] == 2) { 63 //One true neighbour 64 65 //Get index of the one neighbour 66 i=0; k0 = k; 67 while (i<3 && k0==k) { 68 k0 = surrogate_neighbours[index3 + i]; 69 i++; 70 } 71 if (k0 == k) return -1; 72 73 k1 = k; //self 74 75 //Get data 76 q0 = centroid_values[k0]; 77 q1 = centroid_values[k1]; 78 79 x0 = centroids[k0*2]; y0 = centroids[k0*2+1]; 80 x1 = centroids[k1*2]; y1 = centroids[k1*2+1]; 81 82 //Two point gradient 83 _gradient2(x0, y0, x1, y1, q0, q1, &a[k], &b[k]); 84 85 86 //Old (wrong code) 87 //det = x0*y1 - x1*y0; 88 //if (det != 0.0) { 89 // a[k] = (y1*q0 - y0*q1)/det; 90 // b[k] = (x0*q1 - x1*q0)/det; 91 //} 92 } 93 // else: 94 // #No true neighbours - 95 // #Fall back to first order scheme 96 // pass 97 } 98 return 0; 99 99 } 100 100 … … 108 108 double* b) { 109 109 110 111 112 113 for (k=0; k<N; k++){114 115 116 117 118 119 120 121 122 123 124 125 126 127 128 129 130 131 132 133 134 135 136 110 int k, k2, k3, k6; 111 double x, y, x0, y0, x1, y1, x2, y2; 112 113 for (k=0; k<N; k++){ 114 k6 = 6*k; 115 k3 = 3*k; 116 k2 = 2*k; 117 118 //Centroid coordinates 119 x = centroids[k2]; y = centroids[k2+1]; 120 121 //vertex coordinates 122 //x0, y0, x1, y1, x2, y2 = X[k,:] 123 x0 = vertex_coordinates[k6 + 0]; 124 y0 = vertex_coordinates[k6 + 1]; 125 x1 = vertex_coordinates[k6 + 2]; 126 y1 = vertex_coordinates[k6 + 3]; 127 x2 = vertex_coordinates[k6 + 4]; 128 y2 = vertex_coordinates[k6 + 5]; 129 130 //Extrapolate 131 vertex_values[k3+0] = centroid_values[k] + a[k]*(x0-x) + b[k]*(y0-y); 132 vertex_values[k3+1] = centroid_values[k] + a[k]*(x1-x) + b[k]*(y1-y); 133 vertex_values[k3+2] = centroid_values[k] + a[k]*(x2-x) + b[k]*(y2-y); 134 135 } 136 return 0; 137 137 } 138 138 … … 144 144 double* edge_values) { 145 145 146 147 148 149 150 151 152 153 154 155 156 157 158 159 160 161 162 146 int k, k3; 147 double q0, q1, q2; 148 149 150 for (k=0; k<N; k++) { 151 k3 = 3*k; 152 153 q0 = vertex_values[k3 + 0]; 154 q1 = vertex_values[k3 + 1]; 155 q2 = vertex_values[k3 + 2]; 156 157 //printf("%f, %f, %f\n", q0, q1, q2); 158 edge_values[k3 + 0] = 0.5*(q1+q2); 159 edge_values[k3 + 1] = 0.5*(q0+q2); 160 edge_values[k3 + 2] = 0.5*(q0+q1); 161 } 162 return 0; 163 163 } 164 164 … … 168 168 double* explicit_update, 169 169 double* semi_implicit_update) { 170 171 172 173 174 175 176 177 178 179 180 181 182 183 184 185 186 187 188 189 190 191 192 193 194 195 196 197 198 199 200 /* for (k=0; k<N; k++) {*/201 /* centroid_values[k] = exp(timestep*semi_implicit_update[k])*centroid_values[k];*/202 /* }*/203 204 205 206 207 208 209 210 211 212 213 214 215 216 170 //Update centroid values based on values stored in 171 //explicit_update and semi_implicit_update as well as given timestep 172 173 174 int k; 175 double denominator, x; 176 177 178 //Divide semi_implicit update by conserved quantity 179 for (k=0; k<N; k++) { 180 x = centroid_values[k]; 181 if (x == 0.0) { 182 semi_implicit_update[k] = 0.0; 183 } else { 184 semi_implicit_update[k] /= x; 185 } 186 } 187 188 189 //Semi implicit updates 190 for (k=0; k<N; k++) { 191 denominator = 1.0 - timestep*semi_implicit_update[k]; 192 if (denominator == 0.0) { 193 return -1; 194 } else { 195 //Update conserved_quantities from semi implicit updates 196 centroid_values[k] /= denominator; 197 } 198 } 199 200 /* for (k=0; k<N; k++) {*/ 201 /* centroid_values[k] = exp(timestep*semi_implicit_update[k])*centroid_values[k];*/ 202 /* }*/ 203 204 205 //Explicit updates 206 for (k=0; k<N; k++) { 207 centroid_values[k] += timestep*explicit_update[k]; 208 } 209 210 211 //MH080605 set semi_impliit_update[k] to 0.0 here, rather than in update_conserved_quantities.py 212 for (k=0;k<N;k++){ 213 semi_implicit_update[k]=0.0; 214 } 215 216 return 0; 217 217 } 218 218 … … 222 222 PyObject *update(PyObject *self, PyObject *args) { 223 223 224 225 226 227 228 229 230 231 232 233 234 235 236 237 238 239 240 241 242 (double*) centroid_values -> data,243 (double*) explicit_update -> data,244 (double*) semi_implicit_update -> data);245 246 247 248 249 250 251 252 253 254 255 256 257 258 224 PyObject *quantity; 225 PyArrayObject *centroid_values, *explicit_update, *semi_implicit_update; 226 227 double timestep; 228 int N, err; 229 230 231 // Convert Python arguments to C 232 if (!PyArg_ParseTuple(args, "Od", &quantity, ×tep)) 233 return NULL; 234 235 centroid_values = get_consecutive_array(quantity, "centroid_values"); 236 explicit_update = get_consecutive_array(quantity, "explicit_update"); 237 semi_implicit_update = get_consecutive_array(quantity, "semi_implicit_update"); 238 239 N = centroid_values -> dimensions[0]; 240 241 err = _update(N, timestep, 242 (double*) centroid_values -> data, 243 (double*) explicit_update -> data, 244 (double*) semi_implicit_update -> data); 245 246 247 if (err != 0) { 248 PyErr_SetString(PyExc_RuntimeError, 249 "Zero division in semi implicit update - call Stephen :)"); 250 return NULL; 251 } 252 253 //Release and return 254 Py_DECREF(centroid_values); 255 Py_DECREF(explicit_update); 256 Py_DECREF(semi_implicit_update); 257 258 return Py_BuildValue(""); 259 259 } 260 260 … … 262 262 PyObject *interpolate_from_vertices_to_edges(PyObject *self, PyObject *args) { 263 263 264 265 266 267 268 269 270 271 272 273 274 275 276 277 278 264 PyObject *quantity; 265 PyArrayObject *vertex_values, *edge_values; 266 267 int N, err; 268 269 // Convert Python arguments to C 270 if (!PyArg_ParseTuple(args, "O", &quantity)) 271 return NULL; 272 273 vertex_values = get_consecutive_array(quantity, "vertex_values"); 274 edge_values = get_consecutive_array(quantity, "edge_values"); 275 276 N = vertex_values -> dimensions[0]; 277 278 err = _interpolate(N, 279 279 (double*) vertex_values -> data, 280 280 (double*) edge_values -> data); 281 281 282 283 284 285 286 287 288 289 290 291 282 if (err != 0) { 283 PyErr_SetString(PyExc_RuntimeError, "Interpolate could not be computed"); 284 return NULL; 285 } 286 287 //Release and return 288 Py_DECREF(vertex_values); 289 Py_DECREF(edge_values); 290 291 return Py_BuildValue(""); 292 292 } 293 293 … … 295 295 PyObject *compute_gradients(PyObject *self, PyObject *args) { 296 296 297 298 299 300 301 302 303 304 305 306 307 308 309 310 311 312 313 314 315 316 317 318 319 320 321 322 323 324 325 326 327 328 329 330 331 332 333 334 335 336 337 338 339 340 341 342 343 344 345 346 347 348 349 350 351 352 353 354 355 356 357 297 PyObject *quantity, *domain, *R; 298 PyArrayObject 299 *centroids, //Coordinates at centroids 300 *centroid_values, //Values at centroids 301 *number_of_boundaries, //Number of boundaries for each triangle 302 *surrogate_neighbours, //True neighbours or - if one missing - self 303 *a, *b; //Return values 304 305 int dimensions[1], N, err; 306 307 // Convert Python arguments to C 308 if (!PyArg_ParseTuple(args, "O", &quantity)) 309 return NULL; 310 311 domain = PyObject_GetAttrString(quantity, "domain"); 312 if (!domain) 313 return NULL; 314 315 //Get pertinent variables 316 317 centroids = get_consecutive_array(domain, "centroid_coordinates"); 318 centroid_values = get_consecutive_array(quantity, "centroid_values"); 319 surrogate_neighbours = get_consecutive_array(domain, "surrogate_neighbours"); 320 number_of_boundaries = get_consecutive_array(domain, "number_of_boundaries"); 321 322 N = centroid_values -> dimensions[0]; 323 324 //Release 325 Py_DECREF(domain); 326 327 //Allocate space for return vectors a and b (don't DECREF) 328 dimensions[0] = N; 329 a = (PyArrayObject *) PyArray_FromDims(1, dimensions, PyArray_DOUBLE); 330 b = (PyArrayObject *) PyArray_FromDims(1, dimensions, PyArray_DOUBLE); 331 332 333 334 err = _compute_gradients(N, 335 (double*) centroids -> data, 336 (double*) centroid_values -> data, 337 (long*) number_of_boundaries -> data, 338 (long*) surrogate_neighbours -> data, 339 (double*) a -> data, 340 (double*) b -> data); 341 342 if (err != 0) { 343 PyErr_SetString(PyExc_RuntimeError, "Gradient could not be computed"); 344 return NULL; 345 } 346 347 //Release 348 Py_DECREF(centroids); 349 Py_DECREF(centroid_values); 350 Py_DECREF(number_of_boundaries); 351 Py_DECREF(surrogate_neighbours); 352 353 //Build result, release and return 354 R = Py_BuildValue("OO", PyArray_Return(a), PyArray_Return(b)); 355 Py_DECREF(a); 356 Py_DECREF(b); 357 return R; 358 358 } 359 359 … … 362 362 PyObject *extrapolate_second_order(PyObject *self, PyObject *args) { 363 363 364 365 366 367 368 369 370 371 372 373 374 375 376 377 378 379 380 381 382 383 384 385 386 387 388 389 390 391 392 393 394 395 396 397 398 399 400 401 402 403 404 405 406 407 408 409 410 411 412 413 414 415 416 417 418 419 420 421 422 423 424 425 426 427 428 429 430 431 432 433 434 435 436 437 438 439 440 441 442 443 444 445 446 447 448 449 450 451 452 364 PyObject *quantity, *domain; 365 PyArrayObject 366 *centroids, //Coordinates at centroids 367 *centroid_values, //Values at centroids 368 *vertex_coordinates, //Coordinates at vertices 369 *vertex_values, //Values at vertices 370 *number_of_boundaries, //Number of boundaries for each triangle 371 *surrogate_neighbours, //True neighbours or - if one missing - self 372 *a, *b; //Gradients 373 374 //int N, err; 375 int dimensions[1], N, err; 376 //double *a, *b; //Gradients 377 378 // Convert Python arguments to C 379 if (!PyArg_ParseTuple(args, "O", &quantity)) 380 return NULL; 381 382 domain = PyObject_GetAttrString(quantity, "domain"); 383 if (!domain) 384 return NULL; 385 386 //Get pertinent variables 387 centroids = get_consecutive_array(domain, "centroid_coordinates"); 388 centroid_values = get_consecutive_array(quantity, "centroid_values"); 389 surrogate_neighbours = get_consecutive_array(domain, "surrogate_neighbours"); 390 number_of_boundaries = get_consecutive_array(domain, "number_of_boundaries"); 391 vertex_coordinates = get_consecutive_array(domain, "vertex_coordinates"); 392 vertex_values = get_consecutive_array(quantity, "vertex_values"); 393 394 N = centroid_values -> dimensions[0]; 395 396 //Release 397 Py_DECREF(domain); 398 399 //Allocate space for return vectors a and b (don't DECREF) 400 dimensions[0] = N; 401 a = (PyArrayObject *) PyArray_FromDims(1, dimensions, PyArray_DOUBLE); 402 b = (PyArrayObject *) PyArray_FromDims(1, dimensions, PyArray_DOUBLE); 403 404 //FIXME: Odd that I couldn't use normal arrays 405 //Allocate space for return vectors a and b (don't DECREF) 406 //a = (double*) malloc(N * sizeof(double)); 407 //if (!a) return NULL; 408 //b = (double*) malloc(N * sizeof(double)); 409 //if (!b) return NULL; 410 411 412 err = _compute_gradients(N, 413 (double*) centroids -> data, 414 (double*) centroid_values -> data, 415 (long*) number_of_boundaries -> data, 416 (long*) surrogate_neighbours -> data, 417 (double*) a -> data, 418 (double*) b -> data); 419 420 if (err != 0) { 421 PyErr_SetString(PyExc_RuntimeError, "Gradient could not be computed"); 422 return NULL; 423 } 424 425 err = _extrapolate(N, 426 (double*) centroids -> data, 427 (double*) centroid_values -> data, 428 (double*) vertex_coordinates -> data, 429 (double*) vertex_values -> data, 430 (double*) a -> data, 431 (double*) b -> data); 432 433 434 if (err != 0) { 435 PyErr_SetString(PyExc_RuntimeError, 436 "Internal function _extrapolate failed"); 437 return NULL; 438 } 439 440 441 442 //Release 443 Py_DECREF(centroids); 444 Py_DECREF(centroid_values); 445 Py_DECREF(number_of_boundaries); 446 Py_DECREF(surrogate_neighbours); 447 Py_DECREF(vertex_coordinates); 448 Py_DECREF(vertex_values); 449 Py_DECREF(a); 450 Py_DECREF(b); 451 452 return Py_BuildValue(""); 453 453 } 454 454 … … 457 457 PyObject *limit(PyObject *self, PyObject *args) { 458 458 459 PyObject *quantity, *domain, *Tmp; 460 PyArrayObject 461 *qv, //Conserved quantities at vertices 462 *qc, //Conserved quantities at centroids 463 *neighbours; 464 465 int k, i, n, N, k3; 466 double beta_w; //Safety factor 467 double *qmin, *qmax, qn; 468 469 // Convert Python arguments to C 470 if (!PyArg_ParseTuple(args, "O", &quantity)) 471 return NULL; 472 473 domain = PyObject_GetAttrString(quantity, "domain"); 474 if (!domain) 475 return NULL; 476 477 //neighbours = (PyArrayObject*) PyObject_GetAttrString(domain, "neighbours"); 478 neighbours = get_consecutive_array(domain, "neighbours"); 479 480 //Get safety factor beta_w 481 Tmp = PyObject_GetAttrString(domain, "beta_w"); 482 if (!Tmp) 483 return NULL; 484 485 beta_w = PyFloat_AsDouble(Tmp); 486 487 Py_DECREF(Tmp); 488 Py_DECREF(domain); 489 490 491 qc = get_consecutive_array(quantity, "centroid_values"); 492 qv = get_consecutive_array(quantity, "vertex_values"); 493 494 495 N = qc -> dimensions[0]; 496 497 //Find min and max of this and neighbour's centroid values 498 qmin = malloc(N * sizeof(double)); 499 qmax = malloc(N * sizeof(double)); 500 for (k=0; k<N; k++) { 501 k3=k*3; 502 503 qmin[k] = ((double*) qc -> data)[k]; 504 qmax[k] = qmin[k]; 505 506 for (i=0; i<3; i++) { 507 n = ((long*) neighbours -> data)[k3+i]; 508 if (n >= 0) { 509 qn = ((double*) qc -> data)[n]; //Neighbour's centroid value 510 511 qmin[k] = min(qmin[k], qn); 512 qmax[k] = max(qmax[k], qn); 513 } 514 } 515 } 516 517 // Call underlying routine 518 _limit(N, beta_w, (double*) qc -> data, (double*) qv -> data, qmin, qmax); 519 520 free(qmin); 521 free(qmax); 522 return Py_BuildValue(""); 459 PyObject *quantity, *domain, *Tmp; 460 PyArrayObject 461 *qv, //Conserved quantities at vertices 462 *qc, //Conserved quantities at centroids 463 *neighbours; 464 465 int k, i, n, N, k3; 466 double beta_w; //Safety factor 467 double *qmin, *qmax, qn; 468 469 // Convert Python arguments to C 470 if (!PyArg_ParseTuple(args, "O", &quantity)) 471 return NULL; 472 473 domain = PyObject_GetAttrString(quantity, "domain"); 474 if (!domain) 475 return NULL; 476 477 //neighbours = (PyArrayObject*) PyObject_GetAttrString(domain, "neighbours"); 478 neighbours = get_consecutive_array(domain, "neighbours"); 479 480 //Get safety factor beta_w 481 Tmp = PyObject_GetAttrString(domain, "beta_w"); 482 if (!Tmp) 483 return NULL; 484 485 beta_w = PyFloat_AsDouble(Tmp); 486 487 Py_DECREF(Tmp); 488 Py_DECREF(domain); 489 490 491 qc = get_consecutive_array(quantity, "centroid_values"); 492 qv = get_consecutive_array(quantity, "vertex_values"); 493 494 495 N = qc -> dimensions[0]; 496 497 //Find min and max of this and neighbour's centroid values 498 qmin = malloc(N * sizeof(double)); 499 qmax = malloc(N * sizeof(double)); 500 for (k=0; k<N; k++) { 501 k3=k*3; 502 503 qmin[k] = ((double*) qc -> data)[k]; 504 qmax[k] = qmin[k]; 505 506 for (i=0; i<3; i++) { 507 n = ((long*) neighbours -> data)[k3+i]; 508 if (n >= 0) { 509 qn = ((double*) qc -> data)[n]; //Neighbour's centroid value 510 511 qmin[k] = min(qmin[k], qn); 512 qmax[k] = max(qmax[k], qn); 513 } 514 //qmin[k] = max(qmin[k],0.5*((double*) qc -> data)[k]); 515 //qmax[k] = min(qmax[k],2.0*((double*) qc -> data)[k]); 516 } 517 } 518 519 // Call underlying routine 520 _limit(N, beta_w, (double*) qc -> data, (double*) qv -> data, qmin, qmax); 521 522 free(qmin); 523 free(qmax); 524 return Py_BuildValue(""); 523 525 } 524 526 … … 527 529 // Method table for python module 528 530 static struct PyMethodDef MethodTable[] = { 529 530 531 532 533 534 535 536 537 531 {"limit", limit, METH_VARARGS, "Print out"}, 532 {"update", update, METH_VARARGS, "Print out"}, 533 {"compute_gradients", compute_gradients, METH_VARARGS, "Print out"}, 534 {"extrapolate_second_order", extrapolate_second_order, 535 METH_VARARGS, "Print out"}, 536 {"interpolate_from_vertices_to_edges", 537 interpolate_from_vertices_to_edges, 538 METH_VARARGS, "Print out"}, 539 {NULL, NULL, 0, NULL} // sentinel 538 540 }; 539 541 -
anuga_core/source/anuga/abstract_2d_finite_volumes/test_quantity.py
r3514 r3678 941 941 quantity.limit() 942 942 943 943 # limited value for beta_w = 0.9 944 944 assert allclose(quantity.vertex_values[1,:], [2.2, 4.9, 4.9]) 945 # limited values for beta_w = 0.5 946 #assert allclose(quantity.vertex_values[1,:], [3.0, 4.5, 4.5]) 945 947 946 948 -
anuga_core/source/anuga/compile_all.py
r3560 r3678 13 13 14 14 os.chdir('..') 15 os.chdir('shallow_water') 16 execfile('..' + os.sep + 'utilities' + os.sep + 'compile.py') 17 18 os.chdir('..') 15 19 os.chdir('mesh_engine') 16 20 execfile('compile.py') -
anuga_core/source/anuga/config.py
r3642 r3678 39 39 max_smallsteps = 50 #Max number of degenerate steps allowed b4 trying first order 40 40 41 manning = 0. 3 #Manning's friction coefficient41 manning = 0.03 #Manning's friction coefficient 42 42 #g = 9.80665 #Gravity 43 43 g = 9.8 … … 67 67 # 68 68 # 69 #There are separate betas for the w-limiter and the h-limiter 70 # 71 # 72 # 69 #There are separate betas for the w, uh, vh and h limiters 73 70 # 74 71 #Good values are: 75 #beta_w = 0.9 76 #beta_h = 0.2 72 beta_w = 0.9 73 beta_w_dry = 0.9 74 beta_uh = 0.9 75 beta_uh_dry = 0.9 76 beta_vh = 0.9 77 beta_vh_dry = 0.9 78 beta_h = 0.2 79 80 # I think these are better SR but they conflict with the unit tests! 81 # beta_w = 1.0 82 # beta_w_dry = 0.2 83 # beta_uh = 1.0 84 # beta_uh_dry = 0.2 85 # beta_vh = 1.0 86 # beta_vh_dry = 0.2 87 # beta_h = 0.2 77 88 78 89 79 80 beta_w = 0.981 beta_h = 0.282 90 CFL = 1.0 #FIXME (ole): Is this in use yet?? 83 91 #(Steve) yes, change domain.CFL to -
anuga_core/source/anuga/examples/netherlands.py
r3659 r3678 15 15 Transmissive_boundary, Constant_height, Constant_stage 16 16 17 from anuga. mesh_factory import rectangular_cross17 from anuga.abstract_2d_finite_volumes.mesh_factory import rectangular_cross 18 18 from Numeric import array 19 #fromvtk_realtime_visualiser import Visualiser19 from anuga.visualiser.vtk_realtime_visualiser import Visualiser 20 20 21 21 class Weir: … … 83 83 N = 130 #size = 33800 84 84 N = 600 #Size = 720000 85 N = 5085 N = 100 86 86 87 87 … … 96 96 97 97 domain.check_integrity() 98 99 #Setup order and all the beta's for the limiters (these should become defaults 98 100 domain.default_order = 2 99 #domain.beta_h=0 101 domain.beta_w = 1.0 102 domain.beta_w_dry = 0.2 103 domain.beta_uh = 1.0 104 domain.beta_uh_dry = 0.2 105 domain.beta_vh = 1.0 106 domain.beta_vh_dry = 0.2 100 107 101 108 #Output params -
anuga_core/source/anuga/examples/sww_file_visualiser_example.py
r3670 r3678 10 10 # The argument to OfflineVisualiser is the path to a sww file 11 11 o = OfflineVisualiser("../../swollen_viewer/tests/cylinders.sww") 12 #o = OfflineVisualiser("../../../../anuga_validation/analytical solutions/circular_second_order.sww") 12 13 13 14 # Specify the height-based-quantities to render. -
anuga_core/source/anuga/shallow_water/__init__.py
r3659 r3678 8 8 # Make selected classes available directly 9 9 from shallow_water_domain import Domain,\ 10 Transmissive_boundary, Reflective_boundary,\ 11 Dirichlet_boundary, Time_boundary, File_boundary,\ 12 Transmissive_Momentum_Set_Stage_boundary,\ 10 Transmissive_boundary, Reflective_boundary,\ 11 Dirichlet_boundary, Time_boundary, File_boundary,\ 12 Transmissive_Momentum_Set_Stage_boundary,\ 13 Dirichlet_Discharge_boundary,\ 13 14 Constant_stage, Constant_height 14 15 15 16 16 -
anuga_core/source/anuga/shallow_water/data_manager.py
r3664 r3678 724 724 725 725 726 #### N BED national exposure database726 #### NEED national exposure database 727 727 728 728 LAT_TITLE = 'LATITUDE' -
anuga_core/source/anuga/shallow_water/shallow_water_domain.py
r3642 r3678 127 127 128 128 129 from anuga.config import minimum_allowed_height, maximum_allowed_speed, g129 from anuga.config import * 130 130 self.minimum_allowed_height = minimum_allowed_height 131 131 self.maximum_allowed_speed = maximum_allowed_speed 132 132 self.g = g 133 self.beta_w = beta_w 134 self.beta_w_dry = beta_w_dry 135 self.beta_uh = beta_uh 136 self.beta_uh_dry = beta_uh_dry 137 self.beta_vh = beta_vh 138 self.beta_vh_dry = beta_vh_dry 139 self.beta_h = beta_h 133 140 134 141 … … 150 157 self.minimum_storable_height = minimum_storable_height 151 158 self.quantities_to_be_stored = ['stage','xmomentum','ymomentum'] 159 152 160 153 161 … … 301 309 302 310 msg = 'Parameter beta_h must be in the interval [0, 1[' 303 assert 0 <= self.beta_h < 1.0, msg311 assert 0 <= self.beta_h <= 1.0, msg 304 312 msg = 'Parameter beta_w must be in the interval [0, 1[' 305 assert 0 <= self.beta_w < 1.0, msg313 assert 0 <= self.beta_w <= 1.0, msg 306 314 307 315 … … 612 620 Xmom = domain.quantities['xmomentum'] 613 621 Ymom = domain.quantities['ymomentum'] 622 Elevation = domain.quantities['elevation'] 614 623 from shallow_water_ext import extrapolate_second_order_sw 615 extrapolate_second_order_sw(domain,domain.surrogate_neighbours, 624 extrapolate_second_order_sw(domain, 625 domain.surrogate_neighbours, 616 626 domain.number_of_boundaries, 617 627 domain.centroid_coordinates, … … 619 629 Xmom.centroid_values, 620 630 Ymom.centroid_values, 631 Elevation.centroid_values, 621 632 domain.vertex_coordinates, 622 633 Stage.vertex_values, 623 634 Xmom.vertex_values, 624 Ymom.vertex_values) 635 Ymom.vertex_values, 636 Elevation.vertex_values) 625 637 626 638 def compute_fluxes_c(domain): … … 699 711 #all of the conserved quantitie 700 712 701 if (domain. order== 1):713 if (domain._order_ == 1): 702 714 for name in domain.conserved_quantities: 703 715 Q = domain.quantities[name] 704 716 Q.extrapolate_first_order() 705 elif domain. order== 2:717 elif domain._order_ == 2: 706 718 domain.extrapolate_second_order_sw() 707 719 else: … … 711 723 for name in domain.conserved_quantities: 712 724 Q = domain.quantities[name] 713 if domain. order== 1:725 if domain._order_ == 1: 714 726 Q.extrapolate_first_order() 715 elif domain. order== 2:727 elif domain._order_ == 2: 716 728 Q.extrapolate_second_order() 717 729 Q.limit() -
anuga_core/source/anuga/shallow_water/shallow_water_ext.c
r3563 r3678 622 622 *xmom_centroid_values, 623 623 *ymom_centroid_values, 624 *elevation_centroid_values, 624 625 *vertex_coordinates, 625 626 *stage_vertex_values, 626 627 *xmom_vertex_values, 627 *ymom_vertex_values; 628 PyObject *domain, *Tmp; 628 *ymom_vertex_values, 629 *elevation_vertex_values; 630 PyObject *domain, *Tmp_w, *Tmp_w_dry, *Tmp_uh, *Tmp_uh_dry, *Tmp_vh, *Tmp_vh_dry; 629 631 //Local variables 630 632 double a, b;//gradient vector, not stored but used to calculate vertex values from centroids … … 632 634 double x,y,x0,y0,x1,y1,x2,y2,xv0,yv0,xv1,yv1,xv2,yv2;//vertices of the auxiliary triangle 633 635 double dx1,dx2,dy1,dy2,dxv0,dxv1,dxv2,dyv0,dyv1,dyv2,dq0,dq1,dq2,area2; 634 double dqv[3], qmin, qmax, beta_w;//provisional jumps from centroids to v'tices and safety factor re limiting 636 double dqv[3], qmin, qmax, hmin; 637 double hc, h0, h1, h2; 638 double beta_w, beta_w_dry, beta_uh, beta_uh_dry, beta_vh, beta_vh_dry, beta_tmp; 639 //provisional jumps from centroids to v'tices and safety factor re limiting 635 640 //by which these jumps are limited 636 641 // Convert Python arguments to C 637 if (!PyArg_ParseTuple(args, "OOOOOOOOOOO ",642 if (!PyArg_ParseTuple(args, "OOOOOOOOOOOOO", 638 643 &domain, 639 644 &surrogate_neighbours, … … 643 648 &xmom_centroid_values, 644 649 &ymom_centroid_values, 650 &elevation_centroid_values, 645 651 &vertex_coordinates, 646 652 &stage_vertex_values, 647 653 &xmom_vertex_values, 648 &ymom_vertex_values)) { 654 &ymom_vertex_values, 655 &elevation_vertex_values)) { 649 656 PyErr_SetString(PyExc_RuntimeError, "Input arguments failed"); 650 657 return NULL; … … 652 659 653 660 //get the safety factor beta_w, set in the config.py file. This is used in the limiting process 654 Tmp = PyObject_GetAttrString(domain, "beta_w"); 655 if (!Tmp) 656 return NULL; 657 beta_w = PyFloat_AsDouble(Tmp); 658 Py_DECREF(Tmp); 661 Tmp_w = PyObject_GetAttrString(domain, "beta_w"); 662 if (!Tmp_w) 663 return NULL; 664 beta_w = PyFloat_AsDouble(Tmp_w); 665 Py_DECREF(Tmp_w); 666 667 Tmp_w_dry = PyObject_GetAttrString(domain, "beta_w_dry"); 668 if (!Tmp_w_dry) 669 return NULL; 670 beta_w_dry = PyFloat_AsDouble(Tmp_w_dry); 671 Py_DECREF(Tmp_w_dry); 672 673 Tmp_uh = PyObject_GetAttrString(domain, "beta_uh"); 674 if (!Tmp_uh) 675 return NULL; 676 beta_uh = PyFloat_AsDouble(Tmp_uh); 677 Py_DECREF(Tmp_uh); 678 679 Tmp_uh_dry = PyObject_GetAttrString(domain, "beta_uh_dry"); 680 if (!Tmp_uh_dry) 681 return NULL; 682 beta_uh_dry = PyFloat_AsDouble(Tmp_uh_dry); 683 Py_DECREF(Tmp_uh_dry); 684 685 Tmp_vh = PyObject_GetAttrString(domain, "beta_vh"); 686 if (!Tmp_vh) 687 return NULL; 688 beta_vh = PyFloat_AsDouble(Tmp_vh); 689 Py_DECREF(Tmp_vh); 690 691 Tmp_vh_dry = PyObject_GetAttrString(domain, "beta_vh_dry"); 692 if (!Tmp_vh_dry) 693 return NULL; 694 beta_vh_dry = PyFloat_AsDouble(Tmp_vh_dry); 695 Py_DECREF(Tmp_vh_dry); 696 659 697 number_of_elements = stage_centroid_values -> dimensions[0]; 660 698 for (k=0; k<number_of_elements; k++) { … … 710 748 //If the mesh is 'weird' near the boundary, the trianlge might be flat or clockwise: 711 749 if (area2<=0) 712 return NULL; 713 750 return NULL; 751 752 //### Calculate heights of neighbouring cells 753 hc = ((double *)stage_centroid_values->data)[k] - ((double *)elevation_centroid_values->data)[k]; 754 h0 = ((double *)stage_centroid_values->data)[k0] - ((double *)elevation_centroid_values->data)[k0]; 755 h1 = ((double *)stage_centroid_values->data)[k1] - ((double *)elevation_centroid_values->data)[k1]; 756 h2 = ((double *)stage_centroid_values->data)[k2] - ((double *)elevation_centroid_values->data)[k2]; 757 hmin = min(hc,min(h0,min(h1,h2))); 714 758 //### stage ### 715 759 //calculate the difference between vertex 0 of the auxiliary triangle and the FV triangle centroid … … 730 774 //and compute jumps from the centroid to the min and max 731 775 find_qmin_and_qmax(dq0,dq1,dq2,&qmin,&qmax); 732 limit_gradient(dqv,qmin,qmax,beta_w);//the gradient will be limited 776 // Playing with dry wet interface 777 hmin = qmin; 778 beta_tmp = beta_w; 779 if (hmin<0.001) 780 beta_tmp = beta_w_dry; 781 limit_gradient(dqv,qmin,qmax,beta_tmp);//the gradient will be limited 733 782 for (i=0;i<3;i++) 734 783 ((double *)stage_vertex_values->data)[k3+i]=((double *)stage_centroid_values->data)[k]+dqv[i]; … … 752 801 //and compute jumps from the centroid to the min and max 753 802 find_qmin_and_qmax(dq0,dq1,dq2,&qmin,&qmax); 754 limit_gradient(dqv,qmin,qmax,beta_w);//the gradient will be limited 803 beta_tmp = beta_uh; 804 if (hmin<0.001) 805 beta_tmp = beta_uh_dry; 806 limit_gradient(dqv,qmin,qmax,beta_tmp);//the gradient will be limited 755 807 for (i=0;i<3;i++) 756 808 ((double *)xmom_vertex_values->data)[k3+i]=((double *)xmom_centroid_values->data)[k]+dqv[i]; … … 774 826 //and compute jumps from the centroid to the min and max 775 827 find_qmin_and_qmax(dq0,dq1,dq2,&qmin,&qmax); 776 limit_gradient(dqv,qmin,qmax,beta_w);//the gradient will be limited 828 beta_tmp = beta_vh; 829 if (hmin<0.001) 830 beta_tmp = beta_vh_dry; 831 limit_gradient(dqv,qmin,qmax,beta_tmp);//the gradient will be limited 777 832 for (i=0;i<3;i++) 778 833 ((double *)ymom_vertex_values->data)[k3+i]=((double *)ymom_centroid_values->data)[k]+dqv[i]; -
anuga_core/source/anuga/shallow_water/test_shallow_water_domain.py
r3632 r3678 1188 1188 1189 1189 1190 domain. order= 11190 domain._order_ = 1 1191 1191 domain.distribute_to_vertices_and_edges() 1192 1192 … … 1234 1234 assert not alltrue(alltrue(greater_equal(L,E-epsilon))) 1235 1235 1236 domain. order= 11236 domain._order_ = 1 1237 1237 domain.distribute_to_vertices_and_edges() 1238 1238 … … 1269 1269 1270 1270 #First order 1271 domain. order= 11271 domain._order_ = 1 1272 1272 domain.distribute_to_vertices_and_edges() 1273 1273 assert allclose(L[1], val1) 1274 1274 1275 1275 #Second order 1276 domain. order= 21276 domain._order_ = 2 1277 1277 domain.distribute_to_vertices_and_edges() 1278 1278 assert allclose(L[1], [2.2, 4.9, 4.9]) … … 1307 1307 assert allclose(b[1], 0.0) 1308 1308 1309 domain. order= 11309 domain._order_ = 1 1310 1310 domain.distribute_to_vertices_and_edges() 1311 1311 assert allclose(L[1], 1.77777778) 1312 1312 1313 domain. order= 21313 domain._order_ = 2 1314 1314 domain.distribute_to_vertices_and_edges() 1315 1315 assert allclose(L[1], [0.57777777, 2.37777778, 2.37777778]) … … 1345 1345 assert allclose(b[1], 3.33333333) 1346 1346 1347 domain. order= 11347 domain._order_ = 1 1348 1348 domain.distribute_to_vertices_and_edges() 1349 1349 assert allclose(L[1], 4.9382716) 1350 1350 1351 domain. order= 21351 domain._order_ = 2 1352 1352 domain.distribute_to_vertices_and_edges() 1353 1353 assert allclose(L[1], [1.07160494, 6.46058131, 7.28262855]) … … 1391 1391 1392 1392 #print E 1393 domain. order= 11393 domain._order_ = 1 1394 1394 domain.distribute_to_vertices_and_edges() 1395 1395 ##assert allclose(L[1], [0.19999999, 20.05, 20.05]) 1396 1396 assert allclose(L[1], [0.1, 20.1, 20.1]) 1397 1397 1398 domain. order= 21398 domain._order_ = 2 1399 1399 domain.distribute_to_vertices_and_edges() 1400 1400 assert allclose(L[1], [0.1, 20.1, 20.1]) … … 1436 1436 1437 1437 #print E 1438 domain. order= 11438 domain._order_ = 1 1439 1439 domain.distribute_to_vertices_and_edges() 1440 1440 ##assert allclose(L[1], [4.19999999, 16.07142857, 20.02857143]) 1441 1441 assert allclose(L[1], [4.1, 16.1, 20.1]) 1442 1442 1443 domain. order= 21443 domain._order_ = 2 1444 1444 domain.distribute_to_vertices_and_edges() 1445 1445 assert allclose(L[1], [4.1, 16.1, 20.1]) … … 1488 1488 1489 1489 #print E 1490 domain. order= 21490 domain._order_ = 2 1491 1491 domain.beta_h = 0.0 #Use first order in h-limiter 1492 1492 domain.distribute_to_vertices_and_edges() … … 2232 2232 domain.smooth = False 2233 2233 domain.visualise = False 2234 domain.default_order=domain. order=22234 domain.default_order=domain._order_=2 2235 2235 2236 2236 # Boundary conditions -
anuga_core/source/anuga/test_all.py
r3573 r3678 76 76 for file in exclude_files: 77 77 print 'WARNING: File '+ file + ' to be excluded from testing' 78 79 80 81 82 83 84 78 try: 79 files.remove(file) 80 except ValueError, e: 81 msg = 'File "%s" was not found in test suite.\n' %file 82 msg += 'Original error is "%s"\n' %e 83 msg += 'Perhaps it should be removed from exclude list?' 84 raise Exception, msg 85 85 86 86 filenameToModuleName = lambda f: os.path.splitext(f)[0] … … 101 101 runner = unittest.TextTestRunner() #verbosity=2 102 102 runner.run(suite) 103 -
anuga_core/source/anuga_parallel/run_parallel_sw_rectangle.py
r3579 r3678 35 35 processor_name = pypar.Get_processor_name() 36 36 37 M = 2237 M = 50 38 38 N = M*numprocs 39 39 … … 69 69 70 70 71 domain.default_order = 1 71 72 72 73 73 74 #Boundaries … … 85 86 """ 86 87 87 def __init__(self, x0=0.25, x1=0.75, y0=0.0, y1=1.0, h=5.0 ):88 def __init__(self, x0=0.25, x1=0.75, y0=0.0, y1=1.0, h=5.0, h0=0.0): 88 89 self.x0 = x0 89 90 self.x1 = x1 … … 91 92 self.y1 = y1 92 93 self.h = h 94 self.h0 = h0 93 95 94 96 def __call__(self, x, y): 95 return self.h *((x>self.x0)&(x<self.x1)&(y>self.y0)&(y<self.y1))97 return self.h0 + self.h*((x>self.x0)&(x<self.x1)&(y>self.y0)&(y<self.y1)) 96 98 97 domain.set_quantity('stage', Set_Stage(0.2, 0.4,0.25, 0.75, 1.0))99 domain.set_quantity('stage', Set_Stage(0.2, 0.4, 0.25, 0.75, 1.0, 0.00)) 98 100 99 101 if myid == 0: … … 106 108 rect = [0.0, 0.0, 1.0, 1.0] 107 109 domain.initialise_visualiser() 110 111 domain.default_order = 2 112 domain.beta_w = 1.0 113 domain.beta_w_dry = 0.2 114 domain.beta_uh = 1.0 115 domain.beta_uh_dry = 0.2 116 domain.beta_vh = 1.0 117 domain.beta_vh_dry = 0.2 108 118 109 119 yieldstep = 0.005 -
anuga_validation/analytical solutions/Analytical_solution_circular_hydraulic_jump.py
r3514 r3678 8 8 9 9 #------------------------------- 10 # Set up path and module imports 11 import sys 12 from os import sep 13 sys.path.append('..'+sep+'pyvolution') 10 # Setup modules 14 11 15 from shallow_water import Domain, Dirichlet_Discharge_boundary16 from shallow_water import Transmissive_Momentum_Set_Stage_boundary, Dirichlet_boundary12 from anuga.shallow_water import Domain, Dirichlet_Discharge_boundary 13 from anuga.shallow_water import Transmissive_Momentum_Set_Stage_boundary, Dirichlet_boundary 17 14 from math import pi, sqrt 18 from mesh_factory import strang_mesh15 from anuga.abstract_2d_finite_volumes.mesh_factory import strang_mesh 19 16 20 17 … … 66 63 #------------------------------------------ 67 64 # Reduction operation for get_vertex_values 68 from anuga. pyvolution.utilimport mean69 domain.reduction = mean65 from anuga.utilities.numerical_tools import mean 66 #domain.reduction = mean 70 67 #domain.reduction = min #Looks better near steep slopes 71 68 … … 116 113 #------------------ 117 114 # Order of accuracy 118 domain.default_order = 1 119 domain.CFL = 0.75 120 #domain.beta_w = 0.5 121 #domain.beta_h = 0.2 122 domain.smooth = True 115 domain.default_order = 2 116 domain.beta_w = 1.0 117 domain.beta_w_dry = 0.2 118 domain.beta_uh = 1.0 119 domain.beta_uh_dry = 0.2 120 domain.beta_vh = 1.0 121 domain.beta_vh_dry = 0.2 122 domain.CFL = 0.5 123 124 #domain.smooth = True 123 125 124 126 … … 135 137 # 136 138 # 139 140 domain.initialise_visualiser() 141 #from anuga.visualiser.vtk_realtime_visualiser import Visualiser 142 137 143 #from realtime_visualisation_new import Visualiser 144 #vis = Visualiser(domain,title="stage") 145 #vis.setup['elevation'] = True 146 #vis.updating['stage'] = True 147 #vis.qcolor['stage'] = (0.0,0.0,0.8) 148 #vis.coloring['stage']= True 138 149 ##vxmom = Visualiser(domain,title='xmomentum',scale_z=10.0) 139 150 ##vymom = Visualiser(domain,title='ymomentum',scale_z=10.0) … … 148 159 149 160 t0 = time.time() 150 for t in domain.evolve(yieldstep = . 1, finaltime = 3.0):161 for t in domain.evolve(yieldstep = .02, finaltime = 3.0): 151 162 domain.write_time() 152 163 #vis.update() 153 164 exp = '(xmomentum**2 + ymomentum**2)**0.5' 154 165 radial_momentum = domain.create_quantity_from_expression(exp) … … 167 178 168 179 # f.write('time = %25.15e wall clock time %g \n' % (domain.time, time.time())) 169 f.write('%10.3f %25.15e %25.15e %25.15e %25.15e \n' % (domain.time, inner_stage, inner_radial_mom, outer_stage, outer_radial_mom))180 # f.write('%10.3f %25.15e %25.15e %25.15e %25.15e \n' % (domain.time, inner_stage, inner_radial_mom, outer_stage, outer_radial_mom)) 170 181 171 182 f.write('time = %25.15e wall clock time %g \n' % (domain.time, time.time()))
Note: See TracChangeset
for help on using the changeset viewer.