Changeset 1007

Timestamp:

Mar 4, 2005, 4:58:03 PM (20 years ago)

Author:

ole

Message:

Added test for balance_deep_and_shallow:
That stage is conserved
Also checked that C-version produces the same as the Python version.
More tests of this would be nice

Location:

inundation/ga/storm_surge/pyvolution

Files:

• shallow_water.py (modified) (1 diff)
• shallow_water_ext.c (modified) (6 diffs)
• test_shallow_water.py (modified) (2 diffs)

inundation/ga/storm_surge/pyvolution/shallow_water.py

@@ -696 +696 @@
     wv = domain.quantities['stage'].vertex_values
     zv = domain.quantities['elevation'].vertex_values
-    hv = wv-zv
+    hv = wv - zv
 
     #Call C-extension
     from shallow_water_ext import h_limiter
-    h_limiter(domain, hc, hv)
-
-    return hv
+    hvbar = h_limiter(domain, hc, hv)
+
+    return hvbar
 
                                             

inundation/ga/storm_surge/pyvolution/shallow_water_ext.c

@@ -239 +239 @@
       alpha = 1.0;  //Flat bed
 
+    //printf("dz = %.3f, alpha = %.8f\n", dz, alpha);
 
     //  Let 
@@ -743 +744 @@
   PyArrayObject 
     *hv, *hc, //Depth at vertices and centroids        
-    //*hvbar,   //Limited depth at vertices (return values) FIXME: Not used?
+    *hvbar,   //Limited depth at vertices (return values) FIXME: Not used?
     *neighbours;
     
   int k, i, n, N, k3;
-  //int dimensions[2];
+  int dimensions[2];
   double beta_h; //Safety factor (see config.py)
   double *hmin, *hmax, hn;
@@ -769 +770 @@
   
   //Create hvbar
-  //dimensions[0] = N;
-  //dimensions[1] = 3;  
-  //hvbar = (PyArrayObject *) PyArray_FromDims(2, dimensions, PyArray_DOUBLE);
+  dimensions[0] = N;
+  dimensions[1] = 3;  
+  hvbar = (PyArrayObject *) PyArray_FromDims(2, dimensions, PyArray_DOUBLE);
  
 
@@ -785 +786 @@
     for (i=0; i<3; i++) {
       n = ((long*) neighbours -> data)[k3+i];
-      //((double*) hvbar -> data)[k3+i] = ((double*) hc -> data)[k];  
+      
+      //Initialise hvbar with valus from hv
+      ((double*) hvbar -> data)[k3+i] = ((double*) hv -> data)[k3+i];  
+      
       if (n >= 0) {
         hn = ((double*) hc -> data)[n]; //Neighbour's centroid value
@@ -796 +800 @@
   
   // Call underlying routine
-  _limit(N, beta_h, (double*) hc -> data, (double*) hv -> data, hmin, hmax);
+  _limit(N, beta_h, (double*) hc -> data, (double*) hvbar -> data, hmin, hmax);
          
   // // //Py_DECREF(domain); //FIXME: NEcessary?         
@@ -803 +807 @@
   
   //return result using PyArray to avoid memory leak  
-  //return PyArray_Return(hvbar);
-  return Py_BuildValue("");    
+  return PyArray_Return(hvbar);
+  //return Py_BuildValue("");    
 }
 

inundation/ga/storm_surge/pyvolution/test_shallow_water.py

@@ -1445 +1445 @@
         assert allclose(X[1,:], [0.0143507718962, 0.0142502147066, 0.00931268339717])
         assert allclose(Y[1,:], [-0.000117062180693, 7.94434448109e-005, -0.000151505429018])
-        
+
+        
+        
+    def test_balance_deep_and_shallow(self):
+        """Test that balanced limiters preserve conserved quantites.
+        """
+        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
+        elements = [ [1,0,2], [1,2,4], [4,2,5], [3,1,4] ]
+
+        mesh = Domain(points, elements)        
+        mesh.check_integrity()
+
+        #Create a deliberate overshoot
+        mesh.set_quantity('stage', [[3,0,3], [2,2,6], [5,3,8], [8,3,5]])
+        mesh.set_quantity('elevation', 0) #Flat bed
+        stage = mesh.quantities['stage']        
+        
+        ref_centroid_values = copy.copy(stage.centroid_values[:]) #Copy
+    
+        #Limit
+        mesh.distribute_to_vertices_and_edges()
+
+        #Assert that quantities are conserved
+        from Numeric import sum
+        for k in range(mesh.number_of_elements):
+            assert allclose (ref_centroid_values[k],
+                             sum(stage.vertex_values[k,:])/3)
+        
+
+        #Now try with a non-flat bed - closely hugging initial stage in places
+        #This will create alphas in the range [0, 0.478260, 1]
+        mesh.set_quantity('stage', [[3,0,3], [2,2,6], [5,3,8], [8,3,5]])
+        mesh.set_quantity('elevation', [[0,0,0], 
+                                        [1.8,1.9,5.9], 
+                                        [4.6,0,0], 
+                                        [0,2,4]])
+        stage = mesh.quantities['stage']        
+        
+        ref_centroid_values = copy.copy(stage.centroid_values[:]) #Copy
+        ref_vertex_values = copy.copy(stage.vertex_values[:]) #Copy     
+    
+        #Limit
+        mesh.distribute_to_vertices_and_edges()
+
+        
+        #Assert that all vertex quantities have changed
+        for k in range(mesh.number_of_elements):
+            #print ref_vertex_values[k,:], stage.vertex_values[k,:]     
+            assert not allclose (ref_vertex_values[k,:], stage.vertex_values[k,:])              
+        #and assert that quantities are still conserved
+        from Numeric import sum
+        for k in range(mesh.number_of_elements):
+            assert allclose (ref_centroid_values[k],
+                             sum(stage.vertex_values[k,:])/3)
+        
+    
+        #Also check that Python and C version produce the same
+        assert allclose (stage.vertex_values,
+                         [[2,2,2],
+                          [1.93333333, 2.03333333, 6.03333333],
+                          [6.93333333, 4.53333333, 4.53333333],
+                          [5.33333333, 5.33333333, 5.33333333]])
+                          
+                         
+
+        
         
     def test_second_order_flat_bed_onestep(self):
@@ -2693 +2769 @@
 if __name__ == "__main__":
     suite = unittest.makeSuite(TestCase,'test')
-    #suite = unittest.makeSuite(TestCase,'test_spatio_temporal_boundary_2')
+    #suite = unittest.makeSuite(TestCase,'test_balance_deep_and_shallow')
     runner = unittest.TextTestRunner()
     runner.run(suite)