[7562] | 1 | #!/usr/bin/env python |
---|
| 2 | # |
---|
| 3 | |
---|
| 4 | # This file was reverted from changeset:5484 to changeset:5470 on 10th July |
---|
| 5 | # by Ole. |
---|
| 6 | |
---|
| 7 | import unittest |
---|
| 8 | import copy |
---|
| 9 | import os |
---|
| 10 | import numpy as num |
---|
[7762] | 11 | |
---|
| 12 | |
---|
| 13 | import sww |
---|
[7562] | 14 | |
---|
| 15 | |
---|
| 16 | class Test_read_sww(unittest.TestCase): |
---|
| 17 | # Class variable |
---|
| 18 | verbose = False |
---|
| 19 | |
---|
| 20 | def set_verbose(self): |
---|
| 21 | Test_read_sww.verbose = True |
---|
| 22 | |
---|
| 23 | def setUp(self): |
---|
| 24 | pass |
---|
| 25 | |
---|
| 26 | def tearDown(self): |
---|
| 27 | pass |
---|
| 28 | |
---|
| 29 | def test_read_sww(self): |
---|
| 30 | """ |
---|
| 31 | Save to an sww file and then read back the info. |
---|
| 32 | Here we store the info "uniquely" |
---|
| 33 | """ |
---|
| 34 | |
---|
| 35 | #--------------------------------------------------------------------- |
---|
| 36 | # Import necessary modules |
---|
| 37 | #--------------------------------------------------------------------- |
---|
[7762] | 38 | from anuga.abstract_2d_finite_volumes.mesh_factory import \ |
---|
| 39 | rectangular_cross |
---|
| 40 | from anuga.shallow_water.shallow_water_domain import Domain |
---|
[7731] | 41 | from boundaries import Reflective_boundary |
---|
| 42 | from anuga.abstract_2d_finite_volumes.generic_boundary_conditions\ |
---|
[7762] | 43 | import Dirichlet_boundary, Time_boundary |
---|
[7562] | 44 | |
---|
| 45 | #--------------------------------------------------------------------- |
---|
| 46 | # Setup computational domain |
---|
| 47 | #--------------------------------------------------------------------- |
---|
| 48 | length = 8. |
---|
| 49 | width = 4. |
---|
| 50 | dx = dy = 2 # Resolution: Length of subdivisions on both axes |
---|
| 51 | |
---|
| 52 | inc = 0.05 # Elevation increment |
---|
| 53 | |
---|
| 54 | points, vertices, boundary = rectangular_cross(int(length/dx), |
---|
| 55 | int(width/dy), |
---|
| 56 | len1=length, |
---|
| 57 | len2=width) |
---|
| 58 | domain = Domain(points, vertices, boundary) |
---|
| 59 | domain.set_name('read_sww_test'+str(domain.processor)) # Output name |
---|
| 60 | domain.set_quantities_to_be_stored({'elevation': 2, |
---|
| 61 | 'stage': 2, |
---|
| 62 | 'xmomentum': 2, |
---|
| 63 | 'ymomentum': 2, |
---|
| 64 | 'friction': 1}) |
---|
| 65 | |
---|
| 66 | domain.set_store_vertices_uniquely(True) |
---|
| 67 | |
---|
| 68 | #--------------------------------------------------------------------- |
---|
| 69 | # Setup initial conditions |
---|
| 70 | #--------------------------------------------------------------------- |
---|
| 71 | domain.set_quantity('elevation', 0.0) # Flat bed initially |
---|
| 72 | domain.set_quantity('friction', 0.01) # Constant friction |
---|
| 73 | domain.set_quantity('stage', 0.0) # Dry initial condition |
---|
| 74 | |
---|
| 75 | #------------------------------------------------------------------ |
---|
| 76 | # Setup boundary conditions |
---|
| 77 | #------------------------------------------------------------------ |
---|
| 78 | Bi = Dirichlet_boundary([0.4, 0, 0]) # Inflow |
---|
| 79 | Br = Reflective_boundary(domain) # Solid reflective wall |
---|
| 80 | Bo = Dirichlet_boundary([-5, 0, 0]) # Outflow |
---|
| 81 | |
---|
| 82 | domain.set_boundary({'left': Bi, 'right': Bo, 'top': Br, 'bottom': Br}) |
---|
| 83 | |
---|
| 84 | #------------------------------------------------------------------- |
---|
| 85 | # Evolve system through time |
---|
| 86 | #------------------------------------------------------------------- |
---|
| 87 | |
---|
| 88 | for t in domain.evolve(yieldstep=1, finaltime=4.0): |
---|
| 89 | pass |
---|
| 90 | |
---|
| 91 | |
---|
| 92 | # Check that quantities have been stored correctly |
---|
| 93 | source = domain.get_name() + '.sww' |
---|
| 94 | |
---|
| 95 | |
---|
| 96 | #x = fid.variables['x'][:] |
---|
| 97 | #y = fid.variables['y'][:] |
---|
| 98 | #stage = fid.variables['stage'][:] |
---|
| 99 | #elevation = fid.variables['elevation'][:] |
---|
| 100 | #fid.close() |
---|
| 101 | |
---|
| 102 | #assert len(stage.shape) == 2 |
---|
| 103 | #assert len(elevation.shape) == 2 |
---|
| 104 | |
---|
| 105 | #M, N = stage.shape |
---|
| 106 | |
---|
[7762] | 107 | sww_file = sww.Read_sww(source) |
---|
[7562] | 108 | |
---|
| 109 | #print 'last frame number',sww_file.get_last_frame_number() |
---|
| 110 | |
---|
| 111 | assert num.allclose(sww_file.x, domain.get_vertex_coordinates()[:,0]) |
---|
| 112 | assert num.allclose(sww_file.y, domain.get_vertex_coordinates()[:,1]) |
---|
| 113 | |
---|
| 114 | |
---|
| 115 | assert num.allclose(sww_file.time, [0.0, 1.0, 2.0, 3.0, 4.0]) |
---|
| 116 | |
---|
| 117 | M = domain.get_number_of_triangles() |
---|
| 118 | |
---|
| 119 | assert num.allclose(num.reshape(num.arange(3*M), (M,3)), sww_file.vertices) |
---|
| 120 | |
---|
| 121 | last_frame_number = sww_file.get_last_frame_number() |
---|
| 122 | assert last_frame_number == 4 |
---|
| 123 | |
---|
| 124 | assert num.allclose(sww_file.get_bounds(), [0.0, length, 0.0, width]) |
---|
| 125 | |
---|
| 126 | assert 'stage' in sww_file.quantities.keys() |
---|
| 127 | assert 'friction' in sww_file.quantities.keys() |
---|
| 128 | assert 'elevation' in sww_file.quantities.keys() |
---|
| 129 | assert 'xmomentum' in sww_file.quantities.keys() |
---|
| 130 | assert 'ymomentum' in sww_file.quantities.keys() |
---|
| 131 | |
---|
| 132 | |
---|
| 133 | for qname, q in sww_file.read_quantities(last_frame_number).items(): |
---|
| 134 | assert num.allclose(domain.get_quantity(qname).get_values(), q) |
---|
| 135 | |
---|
| 136 | #----------------------------------------- |
---|
| 137 | # Start the evolution off again at frame 3 |
---|
| 138 | #----------------------------------------- |
---|
| 139 | sww_file.read_quantities(last_frame_number-1) |
---|
| 140 | |
---|
| 141 | points, vertices, boundary = rectangular_cross(int(length/dx), |
---|
| 142 | int(width/dy), |
---|
| 143 | len1=length, |
---|
| 144 | len2=width) |
---|
| 145 | new_domain = Domain(points, vertices, boundary) |
---|
| 146 | new_domain.set_quantities_to_be_stored(None) |
---|
| 147 | |
---|
| 148 | new_domain.set_store_vertices_uniquely(True) |
---|
| 149 | |
---|
| 150 | for qname, q in sww_file.read_quantities(last_frame_number-1).items(): |
---|
| 151 | new_domain.set_quantity(qname, q) |
---|
| 152 | |
---|
| 153 | |
---|
| 154 | new_domain.set_starttime(sww_file.get_time()) |
---|
| 155 | #------------------------------------------------------------------ |
---|
| 156 | # Setup boundary conditions |
---|
| 157 | #------------------------------------------------------------------ |
---|
| 158 | Bi = Dirichlet_boundary([0.4, 0, 0]) # Inflow |
---|
| 159 | Br = Reflective_boundary(new_domain) # Solid reflective wall |
---|
| 160 | Bo = Dirichlet_boundary([-5, 0, 0]) # Outflow |
---|
| 161 | |
---|
| 162 | new_domain.set_boundary({'left': Bi, 'right': Bo, 'top': Br, 'bottom': Br}) |
---|
| 163 | |
---|
| 164 | #------------------------------------------------------------------- |
---|
| 165 | # Evolve system through time |
---|
| 166 | #------------------------------------------------------------------- |
---|
| 167 | |
---|
| 168 | for t in new_domain.evolve(yieldstep=1.0, finaltime=1.0): |
---|
| 169 | pass |
---|
| 170 | |
---|
| 171 | # Compare new_domain and domain quantities |
---|
| 172 | for quantity in domain.get_quantity_names(): |
---|
| 173 | dv = domain.get_quantity(quantity).get_values() |
---|
| 174 | ndv = new_domain.get_quantity(quantity).get_values() |
---|
| 175 | |
---|
| 176 | assert num.allclose( dv, ndv) |
---|
| 177 | |
---|
| 178 | # Clean up |
---|
| 179 | os.remove(source) |
---|
| 180 | |
---|
| 181 | |
---|
| 182 | if __name__ == "__main__": |
---|
| 183 | suite = unittest.makeSuite(Test_read_sww, 'test') |
---|
| 184 | runner = unittest.TextTestRunner() #verbosity=2) |
---|
| 185 | runner.run(suite) |
---|
| 186 | |
---|