Changeset 221 for inundation/ga/storm_surge/pyvolution/test_shallow_water.py

Timestamp:

Aug 25, 2004, 5:26:37 PM (21 years ago)

Author:

ole

Message:

Painful testing of second order extrapolation.
One bug found: The protection against infinitesimal heighs at centroid. This was not done in the previous version
and was consequently removed.

File:

• inundation/ga/storm_surge/pyvolution/test_shallow_water.py (modified) (1 diff)

inundation/ga/storm_surge/pyvolution/test_shallow_water.py

@@ -960 +960 @@
         #assert allclose(domain.max_timestep, 0.0396825396825) 
         #print domain.quantities['level'].centroid_values
+
+
+    def test_flatbed_first_order(self):
+        from mesh_factory import rectangular
+        from shallow_water import Domain,\
+             Reflective_boundary, Dirichlet_boundary
+
+        from Numeric import array
+
+        #Create basic mesh
+        N = 8
+        points, vertices, boundary = rectangular(N, N)
+
+        #Create shallow water domain
+        domain = Domain(points, vertices, boundary)
+        domain.smooth = False
+        domain.visualise = False
+        domain.default_order=1
+
+        # Boundary conditions
+        Br = Reflective_boundary(domain)
+        Bd = Dirichlet_boundary([0.2,0.,0.])
+        
+        domain.set_boundary({'left': Bd, 'right': Br, 'top': Br, 'bottom': Br})
+        domain.check_integrity()
+
+  
+        #Evolution
+        for t in domain.evolve(yieldstep = 0.02, finaltime = 0.5):
+            pass
+            #domain.write_time()
+
+        assert allclose(domain.min_timestep, 0.0140413643926)
+        assert allclose(domain.max_timestep, 0.0140947355753)
+
+        for i in range(3):
+            assert allclose(domain.quantities['level'].edge_values[:4,i],
+                            [0.10730244,0.12337617,0.11200126,0.12605666])
+
+            assert allclose(domain.quantities['xmomentum'].edge_values[:4,i],
+                            [0.07610894,0.06901572,0.07284461,0.06819712])
+
+            assert allclose(domain.quantities['ymomentum'].edge_values[:4,i],
+                            [-0.0060238, -0.00157404, -0.00309633, -0.0001637])                         
+        
+    def test_flatbed_second_order(self):
+        from mesh_factory import rectangular
+        from shallow_water import Domain,\
+             Reflective_boundary, Dirichlet_boundary
+
+        from Numeric import array
+
+        #Create basic mesh
+        N = 8
+        points, vertices, boundary = rectangular(N, N)
+
+        #Create shallow water domain
+        domain = Domain(points, vertices, boundary)
+        domain.smooth = False
+        domain.visualise = False
+        domain.default_order=2
+
+        # Boundary conditions
+        Br = Reflective_boundary(domain)
+        Bd = Dirichlet_boundary([0.2,0.,0.])
+
+        domain.set_boundary({'left': Bd, 'right': Br, 'top': Br, 'bottom': Br})
+        domain.check_integrity()
+  
+        #Evolution
+        for t in domain.evolve(yieldstep = 0.01, finaltime = 0.03):
+            pass
+
+
+        assert allclose(domain.min_timestep, 0.0210448446782)
+        assert allclose(domain.max_timestep, 0.0210448446782) 
+
+        #print domain.quantities['level'].vertex_values[:4,0]
+        #print domain.quantities['xmomentum'].vertex_values[:4,0]
+        #print domain.quantities['ymomentum'].vertex_values[:4,0]
+        
+        assert allclose(domain.quantities['level'].vertex_values[:4,0],
+                        [0.00101913,0.05352143,0.00104852,0.05354394]) 
+        assert allclose(domain.quantities['xmomentum'].vertex_values[:4,0],
+                        [0.00090581,0.03685719,0.00088303,0.03687313])
+        assert allclose(domain.quantities['ymomentum'].vertex_values[:4,0],
+                        [-0.00139463,0.0006156,-0.00060364,0.00061827])
+
+
+        
+    def test_flatbed_second_order_distribute(self):
+        #Use real data from pyvolution 2
+        #painfully setup and extracted.
+        from mesh_factory import rectangular
+        from shallow_water import Domain,\
+             Reflective_boundary, Dirichlet_boundary
+
+        from Numeric import array
+
+        #Create basic mesh
+        N = 8
+        points, vertices, boundary = rectangular(N, N)
+
+        #Create shallow water domain
+        domain = Domain(points, vertices, boundary)
+        domain.smooth = False
+        domain.visualise = False
+        domain.default_order=domain.order=2
+
+        # Boundary conditions
+        Br = Reflective_boundary(domain)
+        Bd = Dirichlet_boundary([0.2,0.,0.])
+
+        domain.set_boundary({'left': Bd, 'right': Br, 'top': Br, 'bottom': Br})
+        domain.check_integrity()
+
+
+
+        for V in [False, True]:
+            if V:
+                #Set centroids as if system had been evolved            
+                L = zeros(2*N*N, Float)
+                L[:32] = [7.21205592e-003, 5.35214298e-002, 1.00910824e-002,
+                          5.35439433e-002, 1.00910824e-002, 5.35439433e-002,
+                          1.00910824e-002, 5.35439433e-002, 1.00910824e-002,
+                          5.35439433e-002, 1.00910824e-002, 5.35439433e-002,
+                          1.00910824e-002, 5.35393928e-002, 1.02344264e-002,
+                          5.59605058e-002, 0.00000000e+000, 3.31027800e-004,
+                          0.00000000e+000, 4.37962142e-005, 0.00000000e+000,
+                          4.37962142e-005, 0.00000000e+000, 4.37962142e-005,
+                          0.00000000e+000, 4.37962142e-005, 0.00000000e+000,
+                          4.37962142e-005, 0.00000000e+000, 4.37962142e-005,
+                          0.00000000e+000, 5.57305948e-005]
+        
+                X = zeros(2*N*N, Float)
+                X[:32] = [6.48351607e-003, 3.68571894e-002, 8.50733285e-003,
+                          3.68731327e-002, 8.50733285e-003, 3.68731327e-002,
+                          8.50733285e-003, 3.68731327e-002, 8.50733285e-003,
+                          3.68731327e-002, 8.50733285e-003, 3.68731327e-002,
+                          8.50733285e-003, 3.68693861e-002, 8.65220973e-003,
+                          3.85055387e-002, 0.00000000e+000, 2.86060840e-004,
+                          0.00000000e+000, 3.58905503e-005, 0.00000000e+000,
+                          3.58905503e-005, 0.00000000e+000, 3.58905503e-005,
+                          0.00000000e+000, 3.58905503e-005, 0.00000000e+000,
+                          3.58905503e-005, 0.00000000e+000, 3.58905503e-005,
+                          0.00000000e+000, 4.57662812e-005]
+
+                Y = zeros(2*N*N, Float)
+                Y[:32]=[-1.39463104e-003, 6.15600298e-004, -6.03637382e-004,
+                        6.18272251e-004, -6.03637382e-004, 6.18272251e-004,
+                        -6.03637382e-004, 6.18272251e-004, -6.03637382e-004,
+                        6.18272251e-004, -6.03637382e-004, 6.18272251e-004,
+                        -6.03637382e-004, 6.18599320e-004, -6.74622797e-004,
+                        -1.48934756e-004, 0.00000000e+000, -5.35079969e-005,
+                        0.00000000e+000, -2.57264987e-005, 0.00000000e+000,
+                        -2.57264987e-005, 0.00000000e+000, -2.57264987e-005,
+                        0.00000000e+000, -2.57264987e-005, 0.00000000e+000,
+                        -2.57264987e-005, 0.00000000e+000, -2.57264987e-005,
+                        0.00000000e+000, -2.57635178e-005]
+
+        
+                domain.set_quantity('level', L, 'centroids')
+                domain.set_quantity('xmomentum', X, 'centroids')
+                domain.set_quantity('ymomentum', Y, 'centroids')
+        
+                domain.check_integrity()
+            else:    
+                #Evolution
+                for t in domain.evolve(yieldstep = 0.01, finaltime = 0.03):
+                    pass
+                assert allclose(domain.min_timestep, 0.0210448446782)
+                assert allclose(domain.max_timestep, 0.0210448446782)
+        
+
+            #Centroids were correct but not vertices.
+            #Hence the check of distribute below.
+            assert allclose(domain.quantities['level'].centroid_values[:4],
+                            [0.00721206,0.05352143,0.01009108,0.05354394])
+
+            assert allclose(domain.quantities['xmomentum'].centroid_values[:4],
+                            [0.00648352,0.03685719,0.00850733,0.03687313])
+
+            assert allclose(domain.quantities['ymomentum'].centroid_values[:4],
+                            [-0.00139463,0.0006156,-0.00060364,0.00061827])
+
+            #print 'C17=', domain.quantities['xmomentum'].centroid_values[17]
+            #print 'C19=', domain.quantities['xmomentum'].centroid_values[19]
+
+            assert allclose(domain.quantities['xmomentum'].centroid_values[17],
+                            0.00028606084)
+
+            import copy
+            XX = copy.copy(domain.quantities['xmomentum'].centroid_values)
+            assert allclose(domain.quantities['xmomentum'].centroid_values, XX)
+        
+            domain.distribute_to_vertices_and_edges()
+
+            assert allclose(domain.quantities['xmomentum'].centroid_values, XX)
+
+            assert allclose(domain.quantities['xmomentum'].centroid_values[17],
+                            0.00028606084)
+
+
+            assert allclose(domain.quantities['level'].vertex_values[:4,0],
+                            [0.00101913,0.05352143,0.00104852,0.05354394])
+
+            assert allclose(domain.quantities['ymomentum'].vertex_values[:4,0],
+                            [-0.00139463,0.0006156,-0.00060364,0.00061827])
+
+            #This was the culprit. First triangles vertex 0 had an
+            #x-momentum of 0.0064835 instead of 0.00090581 and
+            #third triangle had 0.00850733 instead of 0.00088303
+            #print domain.quantities['xmomentum'].vertex_values[:4,0]
+        
+            assert allclose(domain.quantities['xmomentum'].vertex_values[:4,0],
+                            [0.00090581,0.03685719,0.00088303,0.03687313])        
+
+