Changeset 6703

Timestamp:

Apr 2, 2009, 4:09:50 PM (16 years ago)

Author:

ole

Message:

Nariman's optimisations of compute_fluxes:

Using memset to initialise three arrays and precomputing denominators.
Simple profile improved compute_fluxes from 4.776s to 4.260s and the overall runtime of evolve from 11.513s to 10.985s. This is an overall speedup exceeding 5%.
Yeeehaha

File:

• anuga_core/source/anuga/shallow_water/shallow_water_ext.c (modified) (12 diffs)

anuga_core/source/anuga/shallow_water/shallow_water_ext.c

@@ -18 +18 @@
 #include "math.h"
 #include <stdio.h>
+#include <string.h>
 
 // Shared code snippets
@@ -26 +27 @@
 
 // Computational function for rotation
+// FIXME: Perhaps inline this and profile
 int _rotate(double *q, double n1, double n2) {
   /*Rotate the momentum component q (q[1], q[2])
@@ -164 +166 @@
 // This is used by flux functions
 // Input parameters uh and h may be modified by this function.
+
+// FIXME: Perhaps inline this and profile
 double _compute_speed(double *uh, 
                       double *h, 
@@ -210 +214 @@
   double v_left, v_right;  
   double s_min, s_max, soundspeed_left, soundspeed_right;
-  double denom, z;
+  double denom, inverse_denominator, z;
 
   // Workspace (allocate once, use many)
@@ -229 +233 @@
   _rotate(q_right_rotated, n1, n2);
 
-  z = (z_left+z_right)/2; // Average elevation values. 
-                          // Even though this will nominally allow for discontinuities 
-                          // in the elevation data, there is currently no numerical 
-                          // support for this so results may be strange near jumps in the bed.
+  z = 0.5*(z_left+z_right); // Average elevation values. 
+                            // Even though this will nominally allow for discontinuities 
+                            // in the elevation data, there is currently no numerical 
+                            // support for this so results may be strange near jumps in the bed.
 
   // Compute speeds in x-direction
@@ -280 +284 @@
     *max_speed = 0.0;
   } else {
+    inverse_denominator = 1.0/denom;
     for (i=0; i<3; i++) {
       edgeflux[i] = s_max*flux_left[i] - s_min*flux_right[i];
       edgeflux[i] += s_max*s_min*(q_right_rotated[i]-q_left_rotated[i]);
-      edgeflux[i] /= denom;
+      edgeflux[i] *= inverse_denominator;
     }
 
@@ -1652 +1657 @@
   // Start computation
   call++; // Flag 'id' of flux calculation for this timestep 
-  
+
+    
   // Set explicit_update to zero for all conserved_quantities.
   // This assumes compute_fluxes called before forcing terms
-  for (k=0; k<number_of_elements; k++) {
-    stage_explicit_update[k]=0.0;
-    xmom_explicit_update[k]=0.0;
-    ymom_explicit_update[k]=0.0;  
-  }
+  memset((char*) stage_explicit_update, 0, number_of_elements*sizeof(double));
+  memset((char*) xmom_explicit_update, 0, number_of_elements*sizeof(double));
+  memset((char*) ymom_explicit_update, 0, number_of_elements*sizeof(double));    
+  
 
   // For all triangles
@@ -1666 +1671 @@
     // Loop through neighbours and compute edge flux for each  
     for (i=0; i<3; i++) {
-      ki = k*3+i; // Linear index (triangle k, edge i)
+      ki = k*3+i; // Linear index to edge i of triangle k
       
       if (already_computed_flux[ki] == call)
@@ -1672 +1677 @@
         continue;
         
-        
+      // Get left hand side values from triangle k, edge i 
       ql[0] = stage_edge_values[ki];
       ql[1] = xmom_edge_values[ki];
@@ -1678 +1683 @@
       zl = bed_edge_values[ki];
 
-      // Quantities at neighbour on nearest face
+      // Get right hand side values either from neighbouring triangle
+      // or from boundary array (Quantities at neighbour on nearest face).
       n = neighbours[ki];
       if (n < 0) {
@@ -1689 +1695 @@
         zr = zl; // Extend bed elevation to boundary
       } else {
-        // Neighbour is a real element
+        // Neighbour is a real triangle
         m = neighbour_edges[ki];
         nm = n*3+m; // Linear index (triangle n, edge m)
@@ -1699 +1705 @@
       }
       
+      // Now we have values for this edge - both from left and right side.
       
       if (optimise_dry_cells) {