Changeset 4685 for anuga_core/source/anuga/shallow_water

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

Location:

anuga_core/source/anuga/shallow_water

Files:

• shallow_water_domain.py (modified) (3 diffs)
• shallow_water_ext.c (modified) (9 diffs)
• test_shallow_water_domain.py (modified) (5 diffs)

anuga_core/source/anuga/shallow_water/shallow_water_domain.py

@@ -102 +102 @@
      beta_uh, beta_uh_dry, beta_vh, beta_vh_dry, tight_slope_limiters
 from anuga.config import alpha_balance
+from anuga.config import optimise_dry_cells
 
 
@@ -164 +165 @@
 
         self.tight_slope_limiters = tight_slope_limiters
+        self.optimise_dry_cells = optimise_dry_cells
 
         self.flux_function = flux_function_central
@@ -889 +891 @@
                                          Ymom.explicit_update,
                                          domain.already_computed_flux,
-                                         domain.max_speed)
+                                         domain.max_speed,
+                                         int(domain.optimise_dry_cells))
 
 

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
     

anuga_core/source/anuga/shallow_water/test_shallow_water_domain.py

@@ -229 +229 @@
         domain.compute_fluxes()
         
-        print domain.get_quantity('stage').explicit_update
+        #print domain.get_quantity('stage').explicit_update
         # FIXME (Ole): TODO the general case
         #assert allclose(domain.get_quantity('stage').explicit_update[1], ........??)
@@ -1306 +1306 @@
 
     #####################################################
-    def test_initial_condition(self):
-        """Test that initial condition is output at time == 0
+
+    def test_flux_optimisation(self):
+        """test_flux_optimisation
+        Test that fluxes are correctly computed using
+        dry and still cell exclusions
         """
 
@@ -1344 +1347 @@
         domain.set_boundary({'exterior': Reflective_boundary(domain)})
 
+
+        #  Check that update arrays are initialised to zero 
+        assert allclose(domain.get_quantity('stage').explicit_update, 0)
+        assert allclose(domain.get_quantity('xmomentum').explicit_update, 0)
+        assert allclose(domain.get_quantity('ymomentum').explicit_update, 0)               
+
+
+        # Get true values
+        domain.optimise_dry_cells = False
+        domain.compute_fluxes()
+        stage_ref = copy.copy(domain.get_quantity('stage').explicit_update)
+        xmom_ref = copy.copy(domain.get_quantity('xmomentum').explicit_update)
+        ymom_ref = copy.copy(domain.get_quantity('ymomentum').explicit_update)       
+
+        # Try with flux optimisation
+        domain.optimise_dry_cells = True
+        domain.compute_fluxes()
+
+        assert allclose(stage_ref, domain.get_quantity('stage').explicit_update)
+        assert allclose(xmom_ref, domain.get_quantity('xmomentum').explicit_update)
+        assert allclose(ymom_ref, domain.get_quantity('ymomentum').explicit_update)
+        
+   
+        
+    def test_initial_condition(self):
+        """test_initial_condition
+        Test that initial condition is output at time == 0 and that
+        computed values change as system evolves
+        """
+
+        from anuga.config import g
+        import copy
+
+        a = [0.0, 0.0]
+        b = [0.0, 2.0]
+        c = [2.0, 0.0]
+        d = [0.0, 4.0]
+        e = [2.0, 2.0]
+        f = [4.0, 0.0]
+
+        points = [a, b, c, d, e, f]
+        #bac, bce, ecf, dbe
+        vertices = [ [1,0,2], [1,2,4], [4,2,5], [3,1,4]]
+
+        domain = Domain(points, vertices)
+
+        #Set up for a gradient of (3,0) at mid triangle (bce)
+        def slope(x, y):
+            return 3*x
+
+        h = 0.1
+        def stage(x,y):
+            return slope(x,y)+h
+
+        domain.set_quantity('elevation', slope)
+        domain.set_quantity('stage', stage)
+
+        # Allow slope limiters to work (FIXME (Ole): Shouldn't this be automatic in ANUGA?)     
+        domain.distribute_to_vertices_and_edges()       
+
+        initial_stage = copy.copy(domain.quantities['stage'].vertex_values)
+
+        domain.set_boundary({'exterior': Reflective_boundary(domain)})
+
+        domain.optimise_dry_cells = True
         #Evolution
-        for t in domain.evolve(yieldstep = 1.0, finaltime = 2.0):
+        for t in domain.evolve(yieldstep = 0.5, finaltime = 2.0):
             stage = domain.quantities['stage'].vertex_values
 
@@ -1352 +1420 @@
             else:
                 assert not allclose(stage, initial_stage)
+
 
         os.remove(domain.get_name() + '.sww')
@@ -4799 +4868 @@
         
 if __name__ == "__main__":
-    suite = unittest.makeSuite(Test_Shallow_Water,'test_flux_computation')    
+
+    suite = unittest.makeSuite(Test_Shallow_Water,'test')    
     #suite = unittest.makeSuite(Test_Shallow_Water,'test_tight_slope_limiters')
     #suite = unittest.makeSuite(Test_Shallow_Water,'test_get_maximum_inundation_from_sww')