Changeset 3703
- Timestamp:
- Oct 5, 2006, 5:50:11 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
r3689 r3703 751 751 752 752 #Yield results 753 754 755 756 #FIXME (Ole, 30 April 2006): Do we need this check?757 758 self.time = finaltime753 if finaltime is not None and self.time >= finaltime: 754 755 if self.time > finaltime: 756 #FIXME (Ole, 30 April 2006): Do we need this check? 757 print 'WARNING (domain.py): time overshot finaltime. Contact Ole.Nielsen@ga.gov.au' 758 self.time = finaltime 759 759 760 760 # Yield final time and stop … … 763 763 764 764 765 765 if self.yieldtime >= yieldstep: 766 766 # Yield (intermediate) time and allow inspection of domain 767 767 -
anuga_core/source/anuga/shallow_water/shallow_water_domain.py
r3689 r3703 144 144 self.beta_h = beta_h 145 145 146 146 self.flux_function = flux_function_central 147 #self.flux_function = flux_function_kinetic 148 147 149 self.forcing_terms.append(manning_friction) 148 150 self.forcing_terms.append(gravity) … … 495 497 #################################### 496 498 # Flux computation 497 def flux_function (normal, ql, qr, zl, zr):499 def flux_function_central(normal, ql, qr, zl, zr): 498 500 """Compute fluxes between volumes for the shallow water wave equation 499 501 cast in terms of w = h+z using the 'central scheme' as described in … … 577 579 578 580 return edgeflux, max_speed 581 582 def erfcc(x): 583 584 from math import fabs, exp 585 586 z=fabs(x) 587 t=1.0/(1.0+0.5*z) 588 result=t*exp(-z*z-1.26551223+t*(1.00002368+t*(.37409196+ 589 t*(.09678418+t*(-.18628806+t*(.27886807+t*(-1.13520398+ 590 t*(1.48851587+t*(-.82215223+t*.17087277))))))))) 591 if x < 0.0: 592 result = 2.0-result 593 594 return result 595 596 def flux_function_kinetic(normal, ql, qr, zl, zr): 597 """Compute fluxes between volumes for the shallow water wave equation 598 cast in terms of w = h+z using the 'central scheme' as described in 599 600 Zhang et. al., Advances in Water Resources, 26(6), 2003, 635-647. 601 602 603 Conserved quantities w, uh, vh are stored as elements 0, 1 and 2 604 in the numerical vectors ql an qr. 605 606 Bed elevations zl and zr. 607 """ 608 609 from anuga.config import g, epsilon 610 from math import sqrt 611 from Numeric import array 612 613 #Align momentums with x-axis 614 q_left = rotate(ql, normal, direction = 1) 615 q_right = rotate(qr, normal, direction = 1) 616 617 z = (zl+zr)/2 #Take average of field values 618 619 w_left = q_left[0] #w=h+z 620 h_left = w_left-z 621 uh_left = q_left[1] 622 623 if h_left < epsilon: 624 u_left = 0.0 #Could have been negative 625 h_left = 0.0 626 else: 627 u_left = uh_left/h_left 628 629 630 w_right = q_right[0] #w=h+z 631 h_right = w_right-z 632 uh_right = q_right[1] 633 634 635 if h_right < epsilon: 636 u_right = 0.0 #Could have been negative 637 h_right = 0.0 638 else: 639 u_right = uh_right/h_right 640 641 vh_left = q_left[2] 642 vh_right = q_right[2] 643 644 soundspeed_left = sqrt(g*h_left) 645 soundspeed_right = sqrt(g*h_right) 646 647 #Maximal wave speed 648 s_max = max(u_left+soundspeed_left, u_right+soundspeed_right, 0) 649 650 #Minimal wave speed 651 s_min = min(u_left-soundspeed_left, u_right-soundspeed_right, 0) 652 653 #Flux computation 654 655 F_left = 0.0 656 F_right = 0.0 657 from math import sqrt, pi, exp 658 if h_left > 0.0: 659 F_left = u_left/sqrt(g*h_left) 660 if h_right > 0.0: 661 F_right = u_right/sqrt(g*h_right) 662 663 edgeflux = array([0.0, 0.0, 0.0]) 664 665 edgeflux[0] = h_left*u_left/2.0*erfcc(-F_left) + \ 666 h_left*sqrt(g*h_left)/2.0/sqrt(pi)*exp(-(F_left**2)) + \ 667 h_right*u_right/2.0*erfcc(F_right) - \ 668 h_right*sqrt(g*h_right)/2.0/sqrt(pi)*exp(-(F_right**2)) 669 670 edgeflux[1] = (h_left*u_left**2 + g/2.0*h_left**2)/2.0*erfcc(-F_left) + \ 671 u_left*h_left*sqrt(g*h_left)/2.0/sqrt(pi)*exp(-(F_left**2)) + \ 672 (h_right*u_right**2 + g/2.0*h_right**2)/2.0*erfcc(F_right) - \ 673 u_right*h_right*sqrt(g*h_right)/2.0/sqrt(pi)*exp(-(F_right**2)) 674 675 edgeflux[2] = vh_left*u_left/2.0*erfcc(-F_left) + \ 676 vh_left*sqrt(g*h_left)/2.0/sqrt(pi)*exp(-(F_left**2)) + \ 677 vh_right*u_right/2.0*erfcc(F_right) - \ 678 vh_right*sqrt(g*h_right)/2.0/sqrt(pi)*exp(-(F_right**2)) 679 680 681 edgeflux = rotate(edgeflux, normal, direction=-1) 682 max_speed = max(abs(s_max), abs(s_min)) 683 684 return edgeflux, max_speed 685 579 686 580 687 … … 648 755 649 756 #Flux computation using provided function 650 edgeflux, max_speed = flux_function(normal, ql, qr, zl, zr)757 edgeflux, max_speed = domain.flux_function(normal, ql, qr, zl, zr) 651 758 flux -= edgeflux * domain.edgelengths[k,i] 652 759 … … 712 819 713 820 timestep = float(sys.maxint) 714 from shallow_water_ext import compute_fluxes 715 domain.timestep = compute_fluxes (timestep, domain.epsilon, domain.g,821 from shallow_water_ext import compute_fluxes_ext_central as compute_fluxes_ext 822 domain.timestep = compute_fluxes_ext(timestep, domain.epsilon, domain.g, 716 823 domain.neighbours, 717 824 domain.neighbour_edges, -
anuga_core/source/anuga/shallow_water/shallow_water_ext.c
r3696 r3703 126 126 127 127 // Computational function for flux computation (using stage w=z+h) 128 int flux_function (double *q_left, double *q_right,128 int flux_function_central(double *q_left, double *q_right, 129 129 double z_left, double z_right, 130 130 double n1, double n2, … … 232 232 return 0; 233 233 } 234 235 double erfcc(double x){ 236 double z,t,result; 237 238 z=fabs(x); 239 t=1.0/(1.0+0.5*z); 240 result=t*exp(-z*z-1.26551223+t*(1.00002368+t*(.37409196+ 241 t*(.09678418+t*(-.18628806+t*(.27886807+t*(-1.13520398+ 242 t*(1.48851587+t*(-.82215223+t*.17087277))))))))); 243 if (x < 0.0) result = 2.0-result; 244 245 return result; 246 } 247 248 249 250 // Computational function for flux computation (using stage w=z+h) 251 int flux_function_kinetic(double *q_left, double *q_right, 252 double z_left, double z_right, 253 double n1, double n2, 254 double epsilon, double g, 255 double *edgeflux, double *max_speed) { 256 257 /*Compute fluxes between volumes for the shallow water wave equation 258 cast in terms of the 'stage', w = h+z using 259 the 'central scheme' as described in 260 261 Zhang et. al., Advances in Water Resources, 26(6), 2003, 635-647. 262 */ 263 264 int i; 265 266 double w_left, h_left, uh_left, vh_left, u_left, F_left; 267 double w_right, h_right, uh_right, vh_right, u_right, F_right; 268 double s_min, s_max, soundspeed_left, soundspeed_right; 269 double z; 270 double q_left_copy[3], q_right_copy[3]; 271 272 273 //Copy conserved quantities to protect from modification 274 for (i=0; i<3; i++) { 275 q_left_copy[i] = q_left[i]; 276 q_right_copy[i] = q_right[i]; 277 } 278 279 //Align x- and y-momentum with x-axis 280 _rotate(q_left_copy, n1, n2); 281 _rotate(q_right_copy, n1, n2); 282 283 z = (z_left+z_right)/2; //Take average of field values 284 285 //Compute speeds in x-direction 286 w_left = q_left_copy[0]; // h+z 287 h_left = w_left-z; 288 uh_left = q_left_copy[1]; 289 290 if (h_left < epsilon) { 291 h_left = 0.0; //Could have been negative 292 u_left = 0.0; 293 } else { 294 u_left = uh_left/h_left; 295 } 296 297 w_right = q_right_copy[0]; 298 h_right = w_right-z; 299 uh_right = q_right_copy[1]; 300 301 if (h_right < epsilon) { 302 h_right = 0.0; //Could have been negative 303 u_right = 0.0; 304 } else { 305 u_right = uh_right/h_right; 306 } 307 308 //Momentum in y-direction 309 vh_left = q_left_copy[2]; 310 vh_right = q_right_copy[2]; 311 312 313 //Maximal and minimal wave speeds 314 soundspeed_left = sqrt(g*h_left); 315 soundspeed_right = sqrt(g*h_right); 316 317 s_max = max(u_left+soundspeed_left, u_right+soundspeed_right); 318 if (s_max < 0.0) s_max = 0.0; 319 320 s_min = min(u_left-soundspeed_left, u_right-soundspeed_right); 321 if (s_min > 0.0) s_min = 0.0; 322 323 324 F_left = 0.0; 325 F_right = 0.0; 326 if (h_left > 0.0) F_left = u_left/sqrt(g*h_left); 327 if (h_right > 0.0) F_right = u_right/sqrt(g*h_right); 328 329 for (i=0; i<3; i++) edgeflux[i] = 0.0; 330 *max_speed = 0.0; 331 332 edgeflux[0] = h_left*u_left/2.0*erfcc(-F_left) + \ 333 h_left*sqrt(g*h_left)/2.0/sqrt(pi)*exp(-(F_left*F_left)) + \ 334 h_right*u_right/2.0*erfcc(F_right) - \ 335 h_right*sqrt(g*h_right)/2.0/sqrt(pi)*exp(-(F_right*F_right)); 336 337 edgeflux[1] = (h_left*u_left*u_left + g/2.0*h_left*h_left)/2.0*erfcc(-F_left) + \ 338 u_left*h_left*sqrt(g*h_left)/2.0/sqrt(pi)*exp(-(F_left*F_left)) + \ 339 (h_right*u_right*u_right + g/2.0*h_right*h_right)/2.0*erfcc(F_right) - \ 340 u_right*h_right*sqrt(g*h_right)/2.0/sqrt(pi)*exp(-(F_right*F_right)); 341 342 edgeflux[2] = vh_left*u_left/2.0*erfcc(-F_left) + \ 343 vh_left*sqrt(g*h_left)/2.0/sqrt(pi)*exp(-(F_left*F_left)) + \ 344 vh_right*u_right/2.0*erfcc(F_right) - \ 345 vh_right*sqrt(g*h_right)/2.0/sqrt(pi)*exp(-(F_right*F_right)); 346 347 //Maximal wavespeed 348 *max_speed = max(fabs(s_max), fabs(s_min)); 349 350 //Rotate back 351 _rotate(edgeflux, n1, -n2); 352 353 return 0; 354 } 355 356 357 234 358 235 359 void _manning_friction(double g, double eps, int N, … … 996 1120 } 997 1121 998 PyObject *compute_fluxes (PyObject *self, PyObject *args) {1122 PyObject *compute_fluxes_ext_central(PyObject *self, PyObject *args) { 999 1123 /*Compute all fluxes and the timestep suitable for all volumes 1000 1124 in domain. 1001 1125 1002 Compute total flux for each conserved quantity using "flux_function "1126 Compute total flux for each conserved quantity using "flux_function_central" 1003 1127 1004 1128 Fluxes across each edge are scaled by edgelengths and summed up … … 1134 1258 normal[1] = ((double *) normals -> data)[ki2+1]; 1135 1259 //Edge flux computation 1136 flux_function (ql, qr, zl, zr,1260 flux_function_central(ql, qr, zl, zr, 1137 1261 normal[0], normal[1], 1138 1262 epsilon, g, … … 1174 1298 } 1175 1299 1300 1301 PyObject *compute_fluxes_ext_kinetic(PyObject *self, PyObject *args) { 1302 /*Compute all fluxes and the timestep suitable for all volumes 1303 in domain. 1304 1305 Compute total flux for each conserved quantity using "flux_function_central" 1306 1307 Fluxes across each edge are scaled by edgelengths and summed up 1308 Resulting flux is then scaled by area and stored in 1309 explicit_update for each of the three conserved quantities 1310 stage, xmomentum and ymomentum 1311 1312 The maximal allowable speed computed by the flux_function for each volume 1313 is converted to a timestep that must not be exceeded. The minimum of 1314 those is computed as the next overall timestep. 1315 1316 Python call: 1317 domain.timestep = compute_fluxes(timestep, 1318 domain.epsilon, 1319 domain.g, 1320 domain.neighbours, 1321 domain.neighbour_edges, 1322 domain.normals, 1323 domain.edgelengths, 1324 domain.radii, 1325 domain.areas, 1326 Stage.edge_values, 1327 Xmom.edge_values, 1328 Ymom.edge_values, 1329 Bed.edge_values, 1330 Stage.boundary_values, 1331 Xmom.boundary_values, 1332 Ymom.boundary_values, 1333 Stage.explicit_update, 1334 Xmom.explicit_update, 1335 Ymom.explicit_update, 1336 already_computed_flux) 1337 1338 1339 Post conditions: 1340 domain.explicit_update is reset to computed flux values 1341 domain.timestep is set to the largest step satisfying all volumes. 1342 1343 1344 */ 1345 1346 1347 PyArrayObject *neighbours, *neighbour_edges, 1348 *normals, *edgelengths, *radii, *areas, 1349 *tri_full_flag, 1350 *stage_edge_values, 1351 *xmom_edge_values, 1352 *ymom_edge_values, 1353 *bed_edge_values, 1354 *stage_boundary_values, 1355 *xmom_boundary_values, 1356 *ymom_boundary_values, 1357 *stage_explicit_update, 1358 *xmom_explicit_update, 1359 *ymom_explicit_update, 1360 *already_computed_flux;//tracks whether the flux across an edge has already been computed 1361 1362 1363 //Local variables 1364 double timestep, max_speed, epsilon, g; 1365 double normal[2], ql[3], qr[3], zl, zr; 1366 double edgeflux[3]; //Work arrays for summing up fluxes 1367 1368 int number_of_elements, k, i, m, n; 1369 int ki, nm=0, ki2; //Index shorthands 1370 static long call=1; 1371 1372 1373 // Convert Python arguments to C 1374 if (!PyArg_ParseTuple(args, "dddOOOOOOOOOOOOOOOOOO", 1375 ×tep, 1376 &epsilon, 1377 &g, 1378 &neighbours, 1379 &neighbour_edges, 1380 &normals, 1381 &edgelengths, &radii, &areas, 1382 &tri_full_flag, 1383 &stage_edge_values, 1384 &xmom_edge_values, 1385 &ymom_edge_values, 1386 &bed_edge_values, 1387 &stage_boundary_values, 1388 &xmom_boundary_values, 1389 &ymom_boundary_values, 1390 &stage_explicit_update, 1391 &xmom_explicit_update, 1392 &ymom_explicit_update, 1393 &already_computed_flux)) { 1394 PyErr_SetString(PyExc_RuntimeError, "Input arguments failed"); 1395 return NULL; 1396 } 1397 number_of_elements = stage_edge_values -> dimensions[0]; 1398 call++;//a static local variable to which already_computed_flux is compared 1399 //set explicit_update to zero for all conserved_quantities. 1400 //This assumes compute_fluxes called before forcing terms 1401 for (k=0; k<number_of_elements; k++) { 1402 ((double *) stage_explicit_update -> data)[k]=0.0; 1403 ((double *) xmom_explicit_update -> data)[k]=0.0; 1404 ((double *) ymom_explicit_update -> data)[k]=0.0; 1405 } 1406 //Loop through neighbours and compute edge flux for each 1407 for (k=0; k<number_of_elements; k++) { 1408 for (i=0; i<3; i++) { 1409 ki = k*3+i; 1410 if (((long *) already_computed_flux->data)[ki]==call)//we've already computed the flux across this edge 1411 continue; 1412 ql[0] = ((double *) stage_edge_values -> data)[ki]; 1413 ql[1] = ((double *) xmom_edge_values -> data)[ki]; 1414 ql[2] = ((double *) ymom_edge_values -> data)[ki]; 1415 zl = ((double *) bed_edge_values -> data)[ki]; 1416 1417 //Quantities at neighbour on nearest face 1418 n = ((long *) neighbours -> data)[ki]; 1419 if (n < 0) { 1420 m = -n-1; //Convert negative flag to index 1421 qr[0] = ((double *) stage_boundary_values -> data)[m]; 1422 qr[1] = ((double *) xmom_boundary_values -> data)[m]; 1423 qr[2] = ((double *) ymom_boundary_values -> data)[m]; 1424 zr = zl; //Extend bed elevation to boundary 1425 } else { 1426 m = ((long *) neighbour_edges -> data)[ki]; 1427 nm = n*3+m; 1428 qr[0] = ((double *) stage_edge_values -> data)[nm]; 1429 qr[1] = ((double *) xmom_edge_values -> data)[nm]; 1430 qr[2] = ((double *) ymom_edge_values -> data)[nm]; 1431 zr = ((double *) bed_edge_values -> data)[nm]; 1432 } 1433 // Outward pointing normal vector 1434 // normal = domain.normals[k, 2*i:2*i+2] 1435 ki2 = 2*ki; //k*6 + i*2 1436 normal[0] = ((double *) normals -> data)[ki2]; 1437 normal[1] = ((double *) normals -> data)[ki2+1]; 1438 //Edge flux computation 1439 flux_function_kinetic(ql, qr, zl, zr, 1440 normal[0], normal[1], 1441 epsilon, g, 1442 edgeflux, &max_speed); 1443 //update triangle k 1444 ((long *) already_computed_flux->data)[ki]=call; 1445 ((double *) stage_explicit_update -> data)[k] -= edgeflux[0]*((double *) edgelengths -> data)[ki]; 1446 ((double *) xmom_explicit_update -> data)[k] -= edgeflux[1]*((double *) edgelengths -> data)[ki]; 1447 ((double *) ymom_explicit_update -> data)[k] -= edgeflux[2]*((double *) edgelengths -> data)[ki]; 1448 //update the neighbour n 1449 if (n>=0){ 1450 ((long *) already_computed_flux->data)[nm]=call; 1451 ((double *) stage_explicit_update -> data)[n] += edgeflux[0]*((double *) edgelengths -> data)[nm]; 1452 ((double *) xmom_explicit_update -> data)[n] += edgeflux[1]*((double *) edgelengths -> data)[nm]; 1453 ((double *) ymom_explicit_update -> data)[n] += edgeflux[2]*((double *) edgelengths -> data)[nm]; 1454 } 1455 ///for (j=0; j<3; j++) { 1456 ///flux[j] -= edgeflux[j]*((double *) edgelengths -> data)[ki]; 1457 ///} 1458 //Update timestep 1459 //timestep = min(timestep, domain.radii[k]/max_speed) 1460 //FIXME: SR Add parameter for CFL condition 1461 if ( ((long *) tri_full_flag -> data)[k] == 1) { 1462 if (max_speed > epsilon) { 1463 timestep = min(timestep, ((double *) radii -> data)[k]/max_speed); 1464 //maxspeed in flux_function is calculated as max(|u+a|,|u-a|) 1465 if (n>=0) 1466 timestep = min(timestep, ((double *) radii -> data)[n]/max_speed); 1467 } 1468 } 1469 } // end for i 1470 //Normalise by area and store for when all conserved 1471 //quantities get updated 1472 ((double *) stage_explicit_update -> data)[k] /= ((double *) areas -> data)[k]; 1473 ((double *) xmom_explicit_update -> data)[k] /= ((double *) areas -> data)[k]; 1474 ((double *) ymom_explicit_update -> data)[k] /= ((double *) areas -> data)[k]; 1475 } //end for k 1476 return Py_BuildValue("d", timestep); 1477 } 1478 1176 1479 PyObject *protect(PyObject *self, PyObject *args) { 1177 1480 // … … 1449 1752 {"rotate", (PyCFunction)rotate, METH_VARARGS | METH_KEYWORDS, "Print out"}, 1450 1753 {"extrapolate_second_order_sw", extrapolate_second_order_sw, METH_VARARGS, "Print out"}, 1451 {"compute_fluxes", compute_fluxes, METH_VARARGS, "Print out"}, 1754 {"compute_fluxes_ext_central", compute_fluxes_ext_central, METH_VARARGS, "Print out"}, 1755 {"compute_fluxes_ext_kinetic", compute_fluxes_ext_kinetic, METH_VARARGS, "Print out"}, 1452 1756 {"gravity", gravity, METH_VARARGS, "Print out"}, 1453 1757 {"manning_friction", manning_friction, METH_VARARGS, "Print out"}, -
anuga_core/source/anuga/shallow_water/shallow_water_kinetic_ext.c
r3698 r3703 1103 1103 normal[1] = ((double *) normals -> data)[ki2+1]; 1104 1104 //Edge flux computation 1105 flux_function_ kinetic(ql, qr, zl, zr,1105 flux_function_central(ql, qr, zl, zr, 1106 1106 normal[0], normal[1], 1107 1107 epsilon, g, -
anuga_core/source/anuga/shallow_water/test_shallow_water_domain.py
r3689 r3703 9 9 10 10 from shallow_water_domain import * 11 11 from shallow_water_domain import flux_function_central as flux_function 12 12 13 13 #Variable windfield implemented using functions
Note: See TracChangeset
for help on using the changeset viewer.