Changeset 4726


Ignore:
Timestamp:
Sep 11, 2007, 2:19:57 PM (17 years ago)
Author:
ole
Message:

Separated flux function into Gateway and Computational part.
This makes the code more readable and also improved the speed from
8.60s to 7.97s in the okushiri performance example on nautilus.
That is an 8.6% improvement for the flux computation and a 2% overall
improvement.

File:
1 edited

Legend:

Unmodified
Added
Removed
  • anuga_core/source/anuga/shallow_water/shallow_water_ext.c

    r4714 r4726  
    188188}
    189189
    190 // Computational function for flux computation (using stage w=z+h)
     190// Innermost flux function (using stage w=z+h)
    191191int flux_function_central(double *q_left, double *q_right,
    192192                          double z_left, double z_right,
     
    213213  double s_min, s_max, soundspeed_left, soundspeed_right;
    214214  double denom, z;
     215 
     216  // FIXME (Ole): Try making these static
    215217  double q_left_rotated[3], q_right_rotated[3];
    216218  double flux_right[3], flux_left[3];
    217219
    218   double h0 = H0*H0; //This ensures a good balance when h approaches H0.
    219                      //But evidence suggests that h0 can be as little as
    220                      //epsilon!
     220  double h0 = H0*H0; // This ensures a good balance when h approaches H0.
     221                     // But evidence suggests that h0 can be as little as
     222                     // epsilon!
    221223 
    222224  // Copy conserved quantities to protect from modification
     
    259261
    260262  s_min = min(u_left-soundspeed_left, u_right-soundspeed_right);
    261    if (s_min > 0.0) s_min = 0.0;
     263  if (s_min > 0.0) s_min = 0.0;
    262264
    263265 
     
    272274
    273275 
    274  
    275276  // Flux computation
    276277  denom = s_max-s_min;
    277   if (denom == 0.0) {
     278  if (denom == 0.0) {  // FIXME (Ole): Try using epsilon here
    278279    for (i=0; i<3; i++) edgeflux[i] = 0.0;
    279280    *max_speed = 0.0;
     
    828829*/
    829830
     831
     832
     833// FIXME (Ole): Split this function up into gateway and computational function
     834// as done with the others (most recently flux 11/9/7).
     835// That will make the code faster and more readable.
     836//
    830837PyObject *extrapolate_second_order_sw(PyObject *self, PyObject *args) {
    831838  /*Compute the vertex values based on a linear reconstruction on each triangle
     
    13661373}
    13671374
     1375
     1376// Computational function for flux computation
     1377double _compute_fluxes_central(int number_of_elements,
     1378                               double timestep,
     1379                               double epsilon,
     1380                               double H0,
     1381                               double g,
     1382                               long* neighbours,
     1383                               long* neighbour_edges,
     1384                               double* normals,
     1385                               double* edgelengths,
     1386                               double* radii,
     1387                               double* areas,
     1388                               long* tri_full_flag,
     1389                               double* stage_edge_values,
     1390                               double* xmom_edge_values,
     1391                               double* ymom_edge_values,
     1392                               double* bed_edge_values,
     1393                               double* stage_boundary_values,
     1394                               double* xmom_boundary_values,
     1395                               double* ymom_boundary_values,
     1396                               double* stage_explicit_update,
     1397                               double* xmom_explicit_update,
     1398                               double* ymom_explicit_update,
     1399                               long* already_computed_flux,
     1400                               double* max_speed_array,
     1401                               int optimise_dry_cells) {
     1402                               
     1403  // Local variables
     1404  double max_speed, length, area;
     1405 
     1406  // FIXME (Ole): Try making arrays static 
     1407  double normal[2], ql[3], qr[3], zl, zr;
     1408  double edgeflux[3]; // Work array for summing up fluxes
     1409
     1410  int k, i, m, n;
     1411
     1412  int ki, nm=0, ki2; // Index shorthands
     1413  static long call=1; // Static local variable flagging already computed flux
     1414                               
     1415       
     1416  // Start computation
     1417  call++; // Flag 'id' of flux calculation for this timestep
     1418 
     1419  // Set explicit_update to zero for all conserved_quantities.
     1420  // This assumes compute_fluxes called before forcing terms
     1421  for (k=0; k<number_of_elements; k++) {
     1422    stage_explicit_update[k]=0.0;
     1423    xmom_explicit_update[k]=0.0;
     1424    ymom_explicit_update[k]=0.0; 
     1425  }
     1426
     1427  // For all triangles
     1428  for (k=0; k<number_of_elements; k++) {
     1429   
     1430    // Loop through neighbours and compute edge flux for each 
     1431    for (i=0; i<3; i++) {
     1432      ki = k*3+i; // Linear index (triangle k, edge i)
     1433     
     1434      if (already_computed_flux[ki] == call)
     1435        // We've already computed the flux across this edge
     1436        continue;
     1437       
     1438       
     1439      ql[0] = stage_edge_values[ki];
     1440      ql[1] = xmom_edge_values[ki];
     1441      ql[2] = ymom_edge_values[ki];
     1442      zl =    bed_edge_values[ki];
     1443
     1444      // Quantities at neighbour on nearest face
     1445      n = neighbours[ki];
     1446      if (n < 0) {
     1447        m = -n-1; // Convert negative flag to index
     1448       
     1449        qr[0] = stage_boundary_values[m];
     1450        qr[1] = xmom_boundary_values[m];
     1451        qr[2] = ymom_boundary_values[m];
     1452        zr = zl; // Extend bed elevation to boundary
     1453      } else {
     1454        m = neighbour_edges[ki];
     1455        nm = n*3+m; // Linear index (triangle n, edge m)
     1456       
     1457        qr[0] = stage_edge_values[nm];
     1458        qr[1] = xmom_edge_values[nm];
     1459        qr[2] = ymom_edge_values[nm];
     1460        zr = bed_edge_values[nm];
     1461      }
     1462     
     1463     
     1464      if (optimise_dry_cells) {     
     1465        // Check if flux calculation is necessary across this edge
     1466        // This check will exclude dry cells.
     1467        // This will also optimise cases where zl != zr as
     1468        // long as both are dry
     1469
     1470        if ( fabs(ql[0] - zl) < epsilon &&
     1471             fabs(qr[0] - zr) < epsilon ) {
     1472          // Cell boundary is dry
     1473         
     1474          already_computed_flux[ki] = call; // #k Done 
     1475          if (n>=0)
     1476            already_computed_flux[nm] = call; // #n Done
     1477       
     1478          max_speed = 0.0;
     1479          continue;
     1480        }
     1481      }
     1482   
     1483     
     1484      // Outward pointing normal vector (domain.normals[k, 2*i:2*i+2])
     1485      ki2 = 2*ki; //k*6 + i*2
     1486      normal[0] = normals[ki2];
     1487      normal[1] = normals[ki2+1];
     1488     
     1489
     1490      // Edge flux computation (triangle k, edge i)
     1491      flux_function_central(ql, qr, zl, zr,
     1492                            normal[0], normal[1],
     1493                            epsilon, H0, g,
     1494                            edgeflux, &max_speed);
     1495     
     1496     
     1497      // Multiply edgeflux by edgelength
     1498      length = edgelengths[ki];
     1499      edgeflux[0] *= length;           
     1500      edgeflux[1] *= length;           
     1501      edgeflux[2] *= length;                       
     1502     
     1503     
     1504      // Update triangle k with flux from edge i
     1505      stage_explicit_update[k] -= edgeflux[0];
     1506      xmom_explicit_update[k] -= edgeflux[1];
     1507      ymom_explicit_update[k] -= edgeflux[2];
     1508     
     1509      already_computed_flux[ki] = call; // #k Done
     1510     
     1511     
     1512      // Update neighbour n with same flux but reversed sign
     1513      if (n>=0) {
     1514        stage_explicit_update[n] += edgeflux[0];
     1515        xmom_explicit_update[n] += edgeflux[1];
     1516        ymom_explicit_update[n] += edgeflux[2];
     1517       
     1518        already_computed_flux[nm] = call; // #n Done
     1519      }
     1520
     1521
     1522      // Update timestep based on edge i and possibly neighbour n
     1523      if (tri_full_flag[k] == 1) {
     1524        if (max_speed > epsilon) {
     1525          timestep = min(timestep, radii[k]/max_speed);
     1526          if (n>=0)
     1527            timestep = min(timestep, radii[n]/max_speed);
     1528        }
     1529      }
     1530     
     1531    } // End edge i
     1532   
     1533   
     1534    // Normalise triangle k by area and store for when all conserved
     1535    // quantities get updated
     1536    area = areas[k];
     1537    stage_explicit_update[k] /= area;
     1538    xmom_explicit_update[k] /= area;
     1539    ymom_explicit_update[k] /= area;
     1540   
     1541   
     1542    // Keep track of maximal speeds
     1543    max_speed_array[k] = max_speed;   
     1544   
     1545  } // End triangle k
     1546       
     1547                               
     1548                               
     1549  return timestep;
     1550}
     1551
     1552
    13681553PyObject *compute_fluxes_ext_central(PyObject *self, PyObject *args) {
     1554  /*Compute all fluxes and the timestep suitable for all volumes
     1555    in domain.
     1556
     1557    Compute total flux for each conserved quantity using "flux_function_central"
     1558
     1559    Fluxes across each edge are scaled by edgelengths and summed up
     1560    Resulting flux is then scaled by area and stored in
     1561    explicit_update for each of the three conserved quantities
     1562    stage, xmomentum and ymomentum
     1563
     1564    The maximal allowable speed computed by the flux_function for each volume
     1565    is converted to a timestep that must not be exceeded. The minimum of
     1566    those is computed as the next overall timestep.
     1567
     1568    Python call:
     1569    domain.timestep = compute_fluxes(timestep,
     1570                                     domain.epsilon,
     1571                                     domain.H0,
     1572                                     domain.g,
     1573                                     domain.neighbours,
     1574                                     domain.neighbour_edges,
     1575                                     domain.normals,
     1576                                     domain.edgelengths,
     1577                                     domain.radii,
     1578                                     domain.areas,
     1579                                     tri_full_flag,
     1580                                     Stage.edge_values,
     1581                                     Xmom.edge_values,
     1582                                     Ymom.edge_values,
     1583                                     Bed.edge_values,
     1584                                     Stage.boundary_values,
     1585                                     Xmom.boundary_values,
     1586                                     Ymom.boundary_values,
     1587                                     Stage.explicit_update,
     1588                                     Xmom.explicit_update,
     1589                                     Ymom.explicit_update,
     1590                                     already_computed_flux,
     1591                                     optimise_dry_cells)                                     
     1592
     1593
     1594    Post conditions:
     1595      domain.explicit_update is reset to computed flux values
     1596      domain.timestep is set to the largest step satisfying all volumes.
     1597
     1598
     1599  */
     1600
     1601
     1602  PyArrayObject *neighbours, *neighbour_edges,
     1603    *normals, *edgelengths, *radii, *areas,
     1604    *tri_full_flag,
     1605    *stage_edge_values,
     1606    *xmom_edge_values,
     1607    *ymom_edge_values,
     1608    *bed_edge_values,
     1609    *stage_boundary_values,
     1610    *xmom_boundary_values,
     1611    *ymom_boundary_values,
     1612    *stage_explicit_update,
     1613    *xmom_explicit_update,
     1614    *ymom_explicit_update,
     1615    *already_computed_flux, //Tracks whether the flux across an edge has already been computed
     1616    *max_speed_array; //Keeps track of max speeds for each triangle
     1617
     1618   
     1619  double timestep, epsilon, H0, g;
     1620  int optimise_dry_cells;
     1621   
     1622  // Convert Python arguments to C
     1623  if (!PyArg_ParseTuple(args, "ddddOOOOOOOOOOOOOOOOOOOi",
     1624                        &timestep,
     1625                        &epsilon,
     1626                        &H0,
     1627                        &g,
     1628                        &neighbours,
     1629                        &neighbour_edges,
     1630                        &normals,
     1631                        &edgelengths, &radii, &areas,
     1632                        &tri_full_flag,
     1633                        &stage_edge_values,
     1634                        &xmom_edge_values,
     1635                        &ymom_edge_values,
     1636                        &bed_edge_values,
     1637                        &stage_boundary_values,
     1638                        &xmom_boundary_values,
     1639                        &ymom_boundary_values,
     1640                        &stage_explicit_update,
     1641                        &xmom_explicit_update,
     1642                        &ymom_explicit_update,
     1643                        &already_computed_flux,
     1644                        &max_speed_array,
     1645                        &optimise_dry_cells)) {
     1646    PyErr_SetString(PyExc_RuntimeError, "Input arguments failed");
     1647    return NULL;
     1648  }
     1649
     1650 
     1651  int number_of_elements = stage_edge_values -> dimensions[0];
     1652
     1653  // Call underlying flux computation routine and update
     1654  // the explicit update arrays
     1655  timestep = _compute_fluxes_central(number_of_elements,
     1656                                     timestep,
     1657                                     epsilon,
     1658                                     H0,
     1659                                     g,
     1660                                     (long*) neighbours -> data,                                 
     1661                                     (long*) neighbour_edges -> data,
     1662                                     (double*) normals -> data,
     1663                                     (double*) edgelengths -> data,
     1664                                     (double*) radii  -> data,
     1665                                     (double*) areas  -> data,
     1666                                     (long*) tri_full_flag -> data,
     1667                                     (double*) stage_edge_values -> data,
     1668                                     (double*) xmom_edge_values -> data,
     1669                                     (double*) ymom_edge_values -> data,
     1670                                     (double*) bed_edge_values -> data,
     1671                                     (double*) stage_boundary_values -> data,
     1672                                     (double*) xmom_boundary_values -> data,
     1673                                     (double*) ymom_boundary_values -> data,
     1674                                     (double*) stage_explicit_update -> data,                         
     1675                                     (double*) xmom_explicit_update -> data,
     1676                                     (double*) ymom_explicit_update -> data,
     1677                                     (long*) already_computed_flux -> data,
     1678                                     (double*) max_speed_array -> data,
     1679                                     optimise_dry_cells);
     1680 
     1681  // Return updated flux timestep
     1682  return Py_BuildValue("d", timestep);
     1683}
     1684
     1685
     1686
     1687
     1688// THIS FUNCTION IS NOW OBSOLETE
     1689PyObject *compute_fluxes_ext_central_original(PyObject *self, PyObject *args) {
    13691690  /*Compute all fluxes and the timestep suitable for all volumes
    13701691    in domain.
Note: See TracChangeset for help on using the changeset viewer.