source: trunk/anuga_core/source/anuga/file/test_sww.py @ 7770

Last change on this file since 7770 was 7770, checked in by hudson, 14 years ago

Broke up data_manager into more modules - some failing tests.

File size: 9.1 KB
RevLine 
[7767]1import os
2import unittest
3import tempfile
4import numpy as num
5
6from anuga.coordinate_transforms.geo_reference import Geo_reference
7from csv_file import load_csv_as_array, load_csv_as_dict
8from anuga.abstract_2d_finite_volumes.mesh_factory import rectangular
9from anuga.shallow_water.shallow_water_domain import Domain
[7770]10from sww import load_sww_as_domain, weed, get_mesh_and_quantities_from_file
[7767]11
12# boundary functions
13from anuga.shallow_water.boundaries import Reflective_boundary, \
14            Field_boundary, Transmissive_momentum_set_stage_boundary, \
15            Transmissive_stage_zero_momentum_boundary
16from anuga.abstract_2d_finite_volumes.generic_boundary_conditions\
17     import Transmissive_boundary, Dirichlet_boundary, \
18            Time_boundary, File_boundary, AWI_boundary
19
20
21class Test_sww(unittest.TestCase):
22    def setUp(self):
23        self.verbose = False
24        pass
25
26    def tearDown(self):
27        pass
28       
29    def test_sww2domain1(self):
30        ################################################
31        #Create a test domain, and evolve and save it.
32        ################################################
33        from mesh_factory import rectangular
34
35        #Create basic mesh
36
37        yiel=0.01
38        points, vertices, boundary = rectangular(10,10)
39
40        #Create shallow water domain
41        domain = Domain(points, vertices, boundary)
42        domain.geo_reference = Geo_reference(56,11,11)
43        domain.smooth = False
44        domain.store = True
45        domain.set_name('bedslope')
46        domain.default_order=2
47        #Bed-slope and friction
48        domain.set_quantity('elevation', lambda x,y: -x/3)
49        domain.set_quantity('friction', 0.1)
50        # Boundary conditions
51        from math import sin, pi
52        Br = Reflective_boundary(domain)
53        Bt = Transmissive_boundary(domain)
54        Bd = Dirichlet_boundary([0.2,0.,0.])
55        Bw = Time_boundary(domain=domain,f=lambda t: [(0.1*sin(t*2*pi)), 0.0, 0.0])
56
57        #domain.set_boundary({'left': Bd, 'right': Br, 'top': Br, 'bottom': Br})
58        domain.set_boundary({'left': Bd, 'right': Bd, 'top': Bd, 'bottom': Bd})
59
60        domain.quantities_to_be_stored['xmomentum'] = 2
61        domain.quantities_to_be_stored['ymomentum'] = 2
62        #Initial condition
63        h = 0.05
64        elevation = domain.quantities['elevation'].vertex_values
65        domain.set_quantity('stage', elevation + h)
66
67        domain.check_integrity()
68        #Evolution
69        #domain.tight_slope_limiters = 1
70        for t in domain.evolve(yieldstep = yiel, finaltime = 0.05):
71            #domain.write_time()
72            pass
73
74
75
76        filename = domain.datadir + os.sep + domain.get_name() + '.sww'
77        domain2 = load_sww_as_domain(filename, None, fail_if_NaN=False,
78                                        verbose=self.verbose)
79        #points, vertices, boundary = rectangular(15,15)
80        #domain2.boundary = boundary
81        ###################
82        ##NOW TEST IT!!!
83        ###################
84
85        os.remove(filename)
86
87        bits = ['vertex_coordinates']
88        for quantity in domain.quantities_to_be_stored:
89            bits.append('get_quantity("%s").get_integral()' % quantity)
90            bits.append('get_quantity("%s").get_values()' % quantity)
91
92        for bit in bits:
93            #print 'testing that domain.'+bit+' has been restored'
94            #print bit
95            #print 'done'
96            assert num.allclose(eval('domain.'+bit),eval('domain2.'+bit))
97
98        ######################################
99        #Now evolve them both, just to be sure
100        ######################################x
101        domain.time = 0.
102        from time import sleep
103
104        final = .1
105        domain.set_quantity('friction', 0.1)
106        domain.store = False
107        domain.set_boundary({'left': Bd, 'right': Bd, 'top': Bd, 'bottom': Bd})
108
109
110        for t in domain.evolve(yieldstep = yiel, finaltime = final):
111            #domain.write_time()
112            pass
113
114        final = final - (domain2.starttime-domain.starttime)
115        #BUT since domain1 gets time hacked back to 0:
116        final = final + (domain2.starttime-domain.starttime)
117
118        domain2.smooth = False
119        domain2.store = False
120        domain2.default_order=2
121        domain2.set_quantity('friction', 0.1)
122        #Bed-slope and friction
123        # Boundary conditions
124        Bd2=Dirichlet_boundary([0.2,0.,0.])
125        domain2.boundary = domain.boundary
126        #print 'domain2.boundary'
127        #print domain2.boundary
128        domain2.set_boundary({'left': Bd, 'right': Bd, 'top': Bd, 'bottom': Bd})
129        #domain2.set_boundary({'exterior': Bd})
130
131        domain2.check_integrity()
132
133        for t in domain2.evolve(yieldstep = yiel, finaltime = final):
134            #domain2.write_time()
135            pass
136
137        ###################
138        ##NOW TEST IT!!!
139        ##################
140
141        bits = ['vertex_coordinates']
142
143        for quantity in ['elevation','stage', 'ymomentum','xmomentum']:
144            bits.append('get_quantity("%s").get_integral()' %quantity)
145            bits.append('get_quantity("%s").get_values()' %quantity)
146
147        #print bits
148        for bit in bits:
149            #print bit
150            #print eval('domain.'+bit)
151            #print eval('domain2.'+bit)
152           
153            #print eval('domain.'+bit+'-domain2.'+bit)
154            msg = 'Values in the two domains are different for ' + bit
155            assert num.allclose(eval('domain.'+bit),eval('domain2.'+bit),
156                                rtol=1.e-5, atol=3.e-8), msg
157
158
159
[7770]160    def test_get_mesh_and_quantities_from_sww_file(self):
161        """test_get_mesh_and_quantities_from_sww_file(self):
162        """     
163       
164        # Generate a test sww file with non trivial georeference
165       
166        import time, os
167        from Scientific.IO.NetCDF import NetCDFFile
168
169        # Setup
170        from mesh_factory import rectangular
171
172        # Create basic mesh (100m x 5m)
173        width = 5
174        length = 50
175        t_end = 10
176        points, vertices, boundary = rectangular(length, width, 50, 5)
177
178        # Create shallow water domain
179        domain = Domain(points, vertices, boundary,
180                        geo_reference = Geo_reference(56,308500,6189000))
181
182        domain.set_name('flowtest')
183        swwfile = domain.get_name() + '.sww'
184        domain.set_datadir('.')
185
186        Br = Reflective_boundary(domain)    # Side walls
187        Bd = Dirichlet_boundary([1, 0, 0])  # inflow
188
189        domain.set_boundary( {'left': Bd, 'right': Bd, 'top': Br, 'bottom': Br})
190
191        for t in domain.evolve(yieldstep=1, finaltime = t_end):
192            pass
193
194       
195        # Read it
196
197        # Get mesh and quantities from sww file
198        X = get_mesh_and_quantities_from_file(swwfile,
199                                              quantities=['elevation',
200                                                          'stage',
201                                                          'xmomentum',
202                                                          'ymomentum'], 
203                                              verbose=False)
204        mesh, quantities, time = X
205       
206
207        # Check that mesh has been recovered
208        assert num.alltrue(mesh.triangles == domain.get_triangles())
209        assert num.allclose(mesh.nodes, domain.get_nodes())
210
211        # Check that time has been recovered
212        assert num.allclose(time, range(t_end+1))
213
214        # Check that quantities have been recovered
215        # (sww files use single precision)
216        z=domain.get_quantity('elevation').get_values(location='unique vertices')
217        assert num.allclose(quantities['elevation'], z)
218
219        for q in ['stage', 'xmomentum', 'ymomentum']:
220            # Get quantity at last timestep
221            q_ref=domain.get_quantity(q).get_values(location='unique vertices')
222
223            #print q,quantities[q]
224            q_sww=quantities[q][-1,:]
225
226            msg = 'Quantity %s failed to be recovered' %q
227            assert num.allclose(q_ref, q_sww, atol=1.0e-6), msg
228           
229        # Cleanup
230        os.remove(swwfile)
231       
232       
233
234    def test_weed(self):
235        coordinates1 = [[0.,0.],[1.,0.],[1.,1.],[1.,0.],[2.,0.],[1.,1.]]
236        volumes1 = [[0,1,2],[3,4,5]]
237        boundary1= {(0,1): 'external',(1,2): 'not external',(2,0): 'external',(3,4): 'external',(4,5): 'external',(5,3): 'not external'}
238        coordinates2,volumes2,boundary2=weed(coordinates1,volumes1,boundary1)
239
240        points2 = {(0.,0.):None,(1.,0.):None,(1.,1.):None,(2.,0.):None}
241
242        assert len(points2)==len(coordinates2)
243        for i in range(len(coordinates2)):
244            coordinate = tuple(coordinates2[i])
245            assert points2.has_key(coordinate)
246            points2[coordinate]=i
247
248        for triangle in volumes1:
249            for coordinate in triangle:
250                assert coordinates2[points2[tuple(coordinates1[coordinate])]][0]==coordinates1[coordinate][0]
251                assert coordinates2[points2[tuple(coordinates1[coordinate])]][1]==coordinates1[coordinate][1]
252
[7767]253#################################################################################
254
255if __name__ == "__main__":
256    suite = unittest.makeSuite(Test_sww, 'test')
257    runner = unittest.TextTestRunner(verbosity=1)
258    runner.run(suite)
Note: See TracBrowser for help on using the repository browser.