Changeset 4685 for anuga_core/source/anuga/shallow_water/shallow_water_ext.c

Timestamp:

Aug 28, 2007, 2:23:35 PM (17 years ago)

Author:

ole

Message:

Implemented optimisation excluding dry cells from flux calculation.
All tests pass and the script run_okushiri_profile.py reports an improvement
in compute_fluxes from 11.25s to 8.58s or 24% faster.
The overall computation was about 40s, so this optimisation improved the
total running time for the problem in question by 7%.

This partially addresses ticket:135

File:

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

anuga_core/source/anuga/shallow_water/shallow_water_ext.c

@@ -245 +245 @@
   soundspeed_right = sqrt(g*h_right);
 
-  
   s_max = max(u_left+soundspeed_left, u_right+soundspeed_right);
   if (s_max < 0.0) s_max = 0.0;
 
   s_min = min(u_left-soundspeed_left, u_right-soundspeed_right);
-  if (s_min > 0.0) s_min = 0.0;
-
+   if (s_min > 0.0) s_min = 0.0;
+
+  
   //Flux formulas
   flux_left[0] = u_left*h_left;
@@ -261 +261 @@
   flux_right[2] = u_right*vh_right;
 
-
+  
+  
   //Flux computation
   denom = s_max-s_min;
@@ -1281 +1282 @@
                                      Xmom.explicit_update,
                                      Ymom.explicit_update,
-                                     already_computed_flux)
+                                     already_computed_flux,
+                                     optimise_dry_cells)                                     
 
 
@@ -1308 +1310 @@
     *max_speed_array; //Keeps track of max speeds for each triangle
 
+
   // Local variables
   double timestep, max_speed, epsilon, g, H0, length, area;
+  int optimise_dry_cells=0; // Optimisation flag  
   double normal[2], ql[3], qr[3], zl, zr;
   double edgeflux[3]; // Work array for summing up fluxes
 
-  int number_of_elements, k, i, j, m, n, computation_needed;
+  int number_of_elements, k, i, m, n; //, j, computation_needed;
+
   int ki, nm=0, ki2; // Index shorthands
   static long call=1; // Static local variable flagging already computed flux
@@ -1319 +1324 @@
 
   // Convert Python arguments to C
-  if (!PyArg_ParseTuple(args, "ddddOOOOOOOOOOOOOOOOOOO",
+  if (!PyArg_ParseTuple(args, "ddddOOOOOOOOOOOOOOOOOOOi",
                         &timestep,
                         &epsilon,
@@ -1340 +1345 @@
                         &ymom_explicit_update,
                         &already_computed_flux,
-                        &max_speed_array)) {
+                        &max_speed_array,
+                        &optimise_dry_cells)) {
     PyErr_SetString(PyExc_RuntimeError, "Input arguments failed");
     return NULL;
   }
+  
   
   number_of_elements = stage_edge_values -> dimensions[0];
@@ -1367 +1374 @@
         // We've already computed the flux across this edge
         continue;
+        
         
       ql[0] = ((double *) stage_edge_values -> data)[ki];
@@ -1393 +1401 @@
       
       
-      // Check if flux calculation is necessary across this edge
-      // FIXME (Ole): Work in progress!
-      computation_needed = 0;
-      //for (j=0; j<3; j++) {      
-      //if (ql[j] != qr[j]) computation_needed = 1;
-      //}
-      
-      //if (computation_needed == 0) {
-        //printf("flux exemption identified\n");
+      if (optimise_dry_cells) {      
+        // Check if flux calculation is necessary across this edge
+        // This check will exclude dry cells.
+        // This will also optimise cases where zl != zr as 
+        // long as both are dry
+
+        if ( fabs(ql[0] - zl) < epsilon && 
+             fabs(qr[0] - zr) < epsilon ) {
+          // Cell boundary is dry
+          
+          ((long *) already_computed_flux -> data)[ki] = call; // #k Done       
+          if (n>=0)
+            ((long *) already_computed_flux -> data)[nm] = call; // #n Done
         
-        //((long *) already_computed_flux -> data)[ki] = call; // #k Done       
-        //if (n>=0)
-        //  ((long *) already_computed_flux -> data)[nm] = call; // #n Done
-        
-        //max_speed = 0.0;
-        //continue;
-      //}
+          max_speed = 0.0;
+          continue;
+        }
+      }
       
 
@@ -1459 +1468 @@
         }
       }
+      
     } // End edge i