 Timestamp:
 May 7, 2008, 5:40:11 PM (16 years ago)
 File:

 1 edited
Legend:
 Unmodified
 Added
 Removed

anuga_core/source/anuga/shallow_water/shallow_water_ext.c
r5238 r5290 144 144 145 145 // Compute Froude number; v/sqrt(gh) 146 146 // FIXME (Ole): Not currently in use 147 147 148 double froude_number; 148 149 … … 549 550 double H0, 550 551 int tight_slope_limiters, 552 int use_centroid_velocities, 551 553 double alpha_balance) { 552 554 553 int k, k3, i , excessive_froude_number=0;555 int k, k3, i; //, excessive_froude_number=0; 554 556 555 557 double dz, hmin, alpha, h_diff, hc_k; 556 558 double epsilon = 1.0e6; // FIXME: Temporary measure 557 double g = 9.81; // FIXME: Temporary measure 558 double hv[3], h; // Depths at vertices 559 double Fx, Fy; // Froude numbers 559 double hv[3]; // Depths at vertices 560 //double Fx, Fy; // Froude numbers 560 561 double uc, vc; // Centroid speeds 561 562 … … 701 702 if (alpha < 1) { 702 703 for (i=0; i<3; i++) { 704 703 705 // FIXME (Ole): Simplify when (if) hvbar gets retired 704 706 if (beta_h > epsilon) { … … 709 711 710 712 // Update momentum at vertices 711 712 713 // Update momentum at vertices 714 // Update momentum as a linear combination of 715 // xmomc and ymomc (shallow) and momentum 716 // from extrapolator xmomv and ymomv (deep). 717 718 719 //xmomv[k3+i] = (1alpha)*xmomc[k] + alpha*xmomv[k3+i]; 720 //ymomv[k3+i] = (1alpha)*ymomc[k] + alpha*ymomv[k3+i]; 721 722 if (tight_slope_limiters == 1) { 723 // FIXME(Ole): Here's what I think (as of 17 Nov 2007) 724 // we need to do. Simple and efficient: 713 if (use_centroid_velocities == 1) { 714 // This is a simple, efficient and robust option 715 // It uses first order approximation of velocities, but retains 716 // the order used by stage. 725 717 726 718 // Speeds at centroids … … 738 730 xmomv[k3+i] = uc*hv[i]; 739 731 ymomv[k3+i] = vc*hv[i]; 732 740 733 } else { 741 734 // Update momentum as a linear combination of 742 735 // xmomc and ymomc (shallow) and momentum 743 736 // from extrapolator xmomv and ymomv (deep). 744 // FIXME (Ole): Is this really needed? Could we use the above 745 // instead? 737 // This assumes that values from xmomv and ymomv have 738 // been established e.g. by the gradient limiter. 739 746 740 747 741 xmomv[k3+i] = (1alpha)*xmomc[k] + alpha*xmomv[k3+i]; … … 749 743 750 744 } 751 }752 }753 754 755 756 if (0) { // FIXME(Ole): Disabled while testing balancing of velocities above757 //if (tight_slope_limiters == 1) {758 759 // Ensure that the Froude number is kept realistic at vertices760 // FIXME (Ole): I think it could be used to adjust alpha down761 // whenever Fr is too large. Possible make sure it doesn't deviate762 // too much from the value at the centroid. I like this idea!763 764 // FIXME (Ole): currently only used with tights_SL765 766 // FIXME (Ole): may not be necessary now767 768 excessive_froude_number=0;769 for (i=0; i<3; i++) {770 771 // Recalculate depth at vertex i772 h = wv[k3+i]  zv[k3+i];773 774 Fx = compute_froude_number(xmomv[k3+i], h, g, epsilon);775 Fy = compute_froude_number(ymomv[k3+i], h, g, epsilon);776 777 if ( (fabs(Fx) > 100.0)  (fabs(Fy) > 100.0)) {778 // FIXME: use max_froude  or base it on centroid value of F779 // printf("Excessive Froude number detected: %f or %f\n", Fx, Fy);780 excessive_froude_number=1;781 }782 }783 784 if (excessive_froude_number) {785 786 // printf("Adjusting momentum to first order.\n");787 // Go back to first order (alpha = 0) for this triangle788 for (i=0; i<3; i++) {789 xmomv[k3+i] = xmomc[k];790 ymomv[k3+i] = ymomc[k];791 792 if (beta_h > epsilon) {793 wv[k3+i] = zv[k3+i] + hvbar[k3+i];794 } else {795 wv[k3+i] = zv[k3+i] + hc_k;796 }797 }798 745 } 799 746 } … … 1094 1041 double* ymom_vertex_values, 1095 1042 double* elevation_vertex_values, 1096 int optimise_dry_cells) { 1043 int optimise_dry_cells, 1044 int use_centroid_velocities) { 1097 1045 1098 1046 … … 1359 1307 // One internal neighbour and gradient is in direction of the neighbour's centroid 1360 1308 1361 // Find the only internal neighbour 1309 // Find the only internal neighbour (k1?) 1362 1310 for (k2=k3;k2<k3+3;k2++){ 1363 1311 // Find internal neighbour of triangle k … … 1367 1315 break; 1368 1316 } 1317 1369 1318 if ((k2==k3+3)) { 1370 1319 // If we didn't find an internal neighbour … … 1396 1345 dy2=dx2*dy1; 1397 1346 dx2*=dx1; 1398 1399 1347 1400 1348 … … 1561 1509 double beta_w, beta_w_dry, beta_uh, beta_uh_dry, beta_vh, beta_vh_dry; 1562 1510 double minimum_allowed_height, epsilon; 1563 int optimise_dry_cells, number_of_elements, e ;1511 int optimise_dry_cells, number_of_elements, e, use_centroid_velocities; 1564 1512 1565 1513 // Provisional jumps from centroids to v'tices and safety factor re limiting 1566 1514 // by which these jumps are limited 1567 1515 // Convert Python arguments to C 1568 if (!PyArg_ParseTuple(args, "OOOOOOOOOOOOOi ",1516 if (!PyArg_ParseTuple(args, "OOOOOOOOOOOOOii", 1569 1517 &domain, 1570 1518 &surrogate_neighbours, … … 1580 1528 &ymom_vertex_values, 1581 1529 &elevation_vertex_values, 1582 &optimise_dry_cells)) { 1530 &optimise_dry_cells, 1531 &use_centroid_velocities)) { 1583 1532 1584 1533 PyErr_SetString(PyExc_RuntimeError, … … 1625 1574 (double*) ymom_vertex_values > data, 1626 1575 (double*) elevation_vertex_values > data, 1627 optimise_dry_cells); 1576 optimise_dry_cells, 1577 use_centroid_velocities); 1628 1578 if (e == 1) { 1629 1579 // Use error string set inside computational routine … … 2271 2221 double H0, beta_h; 2272 2222 2273 int N, tight_slope_limiters ; //, err;2223 int N, tight_slope_limiters, use_centroid_velocities; //, err; 2274 2224 2275 2225 // Convert Python arguments to C … … 2318 2268 2319 2269 2270 Tmp = PyObject_GetAttrString(domain, "use_centroid_velocities"); 2271 if (!Tmp) { 2272 PyErr_SetString(PyExc_RuntimeError, "shallow_water_ext.c: balance_deep_and_shallow could not obtain object use_centroid_velocities from domain"); 2273 return NULL; 2274 } 2275 use_centroid_velocities = PyInt_AsLong(Tmp); 2276 Py_DECREF(Tmp); 2277 2278 2279 2320 2280 N = wc > dimensions[0]; 2321 2281 _balance_deep_and_shallow(N, … … 2332 2292 H0, 2333 2293 (int) tight_slope_limiters, 2294 (int) use_centroid_velocities, 2334 2295 alpha_balance); 2335 2296
Note: See TracChangeset
for help on using the changeset viewer.