source: trunk/anuga_core/source/anuga/abstract_2d_finite_volumes/test_generic_boundary_conditions.py @ 8124

Last change on this file since 8124 was 8124, checked in by wilsonr, 13 years ago

Made changes described in ticket 359.

File size: 12.7 KB
Line 
1#!/usr/bin/env python
2"""
3    Generic boundary conditions for a domain.
4   
5    A boundary represents the edge of the model, where inflow, outflow, and
6    reflection can take place.
7   
8    The boundaries in this model can be applied universally across all
9    domain models, without being tied to a particular implementation.
10"""
11
12import unittest
13from math import sqrt, pi
14
15from generic_boundary_conditions import *
16from anuga.abstract_2d_finite_volumes.generic_domain import Generic_Domain
17from anuga.config import epsilon
18
19import numpy as num
20
21
22class Test_Generic_Boundary_Conditions(unittest.TestCase):
23    def setUp(self):
24        pass
25        #print "  Setting up"
26
27    def tearDown(self):
28        pass
29        #print "  Tearing down"
30
31
32    def test_generic(self):
33        b = Boundary()
34
35        try:
36            b.evaluate()
37        except:
38            pass
39        else:
40            raise Exception('Should have raised exception')
41
42
43    def test_dirichlet_empty(self):
44
45        try:
46            Bd = Dirichlet_boundary()
47        except:
48            pass
49        else:
50            raise Exception('Should have raised exception')
51
52    def test_dirichlet(self):
53        x = [3.14,0,0.1]
54        Bd = Dirichlet_boundary(x)
55
56        q = Bd.evaluate()
57        assert num.allclose(q, x)
58
59
60    def test_time(self):
61        from quantity import Quantity
62
63        a = [0.0, 0.0]
64        b = [0.0, 2.0]
65        c = [2.0,0.0]
66        d = [0.0, 4.0]
67        e = [2.0, 2.0]
68        f = [4.0,0.0]
69
70        points = [a, b, c, d, e, f]
71
72        #bac, bce, ecf, dbe
73        elements = [ [1,0,2], [1,2,4], [4,2,5], [3,1,4] ]
74
75        domain = Generic_Domain(points, elements)
76        domain.check_integrity()
77
78        domain.conserved_quantities = ['stage', 'ymomentum']
79        domain.evolved_quantities = ['stage', 'ymomentum']       
80        domain.quantities['stage'] =\
81                                   Quantity(domain, [[1,2,3], [5,5,5],
82                                                     [0,0,9], [-6, 3, 3]])
83
84        domain.quantities['ymomentum'] =\
85                                   Quantity(domain, [[2,3,4], [5,5,5],
86                                                     [0,0,9], [-6, 3, 3]])
87
88
89        domain.check_integrity()
90
91        #Test time bdry, you need to provide a domain and function
92        try:
93            T = Time_boundary(domain)
94        except:
95            pass
96        else:
97            raise Exception('Should have raised exception')
98
99        #Test time bdry, you need to provide a function
100        try:
101            T = Time_boundary()
102        except:
103            pass
104        else:
105            raise Exception('Should have raised exception')
106
107
108        def function(t):
109            return [1.0, 0.0]
110       
111        T = Time_boundary(domain, function)
112
113        from anuga.config import default_boundary_tag
114        domain.set_boundary( {default_boundary_tag: T} )
115
116
117        #FIXME: should not necessarily be true always.
118        #E.g. with None as a boundary object.
119        assert len(domain.boundary) == len(domain.boundary_objects)
120
121        q = T.evaluate(0, 2)  #Vol=0, edge=2
122
123        assert num.allclose(q, [1.0, 0.0])
124
125
126    def test_time_space_boundary(self):
127        from quantity import Quantity
128
129        a = [0.0, 0.0]
130        b = [0.0, 2.0]
131        c = [2.0,0.0]
132        d = [0.0, 4.0]
133        e = [2.0, 2.0]
134        f = [4.0,0.0]
135
136        points = [a, b, c, d, e, f]
137
138        #bac, bce, ecf, dbe
139        elements = [ [1,0,2], [1,2,4], [4,2,5], [3,1,4] ]
140
141        domain = Generic_Domain(points, elements)
142        domain.check_integrity()
143
144        domain.conserved_quantities = ['stage', 'ymomentum']
145        domain.evolved_quantities = ['stage', 'ymomentum']       
146        domain.quantities['stage'] =\
147                                   Quantity(domain, [[1,2,3], [5,5,5],
148                                                     [0,0,9], [-6, 3, 3]])
149
150        domain.quantities['ymomentum'] =\
151                                   Quantity(domain, [[2,3,4], [5,5,5],
152                                                     [0,0,9], [-6, 3, 3]])
153
154
155        domain.check_integrity()
156
157        #Test time space bdry, you need to provide a domain and function
158        try:
159            T = Time_space_boundary(domain)
160        except:
161            pass
162        else:
163            raise Exception('Should have raised exception')
164
165        #Test time bdry, you need to provide a function
166        try:
167            T = Time_space_boundary()
168        except:
169            pass
170        else:
171            raise Exception('Should have raised exception')
172
173
174        def function(t,x,y):
175            return [x,y]
176       
177        T = Time_space_boundary(domain, function)
178
179        from anuga.config import default_boundary_tag
180        domain.set_boundary( {default_boundary_tag: T} )
181
182
183        #FIXME: should not necessarily be true always.
184        #E.g. with None as a boundary object.
185        assert len(domain.boundary) == len(domain.boundary_objects)
186
187        q = T.evaluate(0, 2)  #Vol=0, edge=2
188        assert num.allclose(q, domain.get_edge_midpoint_coordinate(0,2))
189
190
191        q = T.evaluate(1, 1)  #Vol=1, edge=1
192        assert num.allclose(q, domain.get_edge_midpoint_coordinate(1,1))       
193
194
195
196
197
198    def test_transmissive(self):
199        from quantity import Quantity
200
201        a = [0.0, 0.0]
202        b = [0.0, 2.0]
203        c = [2.0,0.0]
204        d = [0.0, 4.0]
205        e = [2.0, 2.0]
206        f = [4.0,0.0]
207
208        points = [a, b, c, d, e, f]
209
210        #bac, bce, ecf, dbe
211        elements = [ [1,0,2], [1,2,4], [4,2,5], [3,1,4] ]
212
213        domain = Generic_Domain(points, elements)
214        domain.check_integrity()
215
216        domain.conserved_quantities = ['stage', 'ymomentum']
217        domain.evolved_quantities = ['stage', 'ymomentum']       
218        domain.quantities['stage'] =\
219                                   Quantity(domain, [[1,2,3], [5,5,5],
220                                                     [0,0,9], [-6, 3, 3]])
221
222        domain.quantities['ymomentum'] =\
223                                   Quantity(domain, [[2,3,4], [5,5,5],
224                                                     [0,0,9], [-6, 3, 3]])
225
226
227        domain.check_integrity()
228
229        #Test transmissve bdry
230        try:
231            T = Transmissive_boundary()
232        except:
233            pass
234        else:
235            raise Exception('Should have raised exception')
236
237        T = Transmissive_boundary(domain)
238
239        from anuga.config import default_boundary_tag
240        domain.set_boundary( {default_boundary_tag: T} )
241
242
243        #FIXME: should not necessarily be true always.
244        #E.g. with None as a boundary object.
245        assert len(domain.boundary) == len(domain.boundary_objects)
246
247        q = T.evaluate(0, 2)  #Vol=0, edge=2
248
249        assert num.allclose(q, [1.5, 2.5])
250
251
252        # Now set the centroid_transmissive_bc flag to true
253        domain.set_centroid_transmissive_bc(True)
254
255        q = T.evaluate(0, 2)  #Vol=0, edge=2
256
257        assert num.allclose(q, [2.0 ,3.0]) # centroid value
258
259
260
261       
262
263    def NOtest_fileboundary_time_only(self):
264        """Test that boundary values can be read from file and interpolated
265        This is using the .tms file format
266       
267        See also test_util for comprenhensive testing of the underlying
268        file_function and also tests in test_datamanager which tests
269        file_function using the sts format
270        """
271        #FIXME (Ole): This test was disabled 18 August 2008 as no
272        # need for this was found. Rather I implemented an Exception
273        # to catch possible errors in the model setup
274       
275        from quantity import Quantity
276        import time, os
277        from math import sin, pi
278        from anuga.config import time_format
279
280        a = [0.0, 0.0]
281        b = [0.0, 2.0]
282        c = [2.0, 0.0]
283        d = [0.0, 4.0]
284        e = [2.0, 2.0]
285        f = [4.0, 0.0]
286
287        points = [a, b, c, d, e, f]
288
289        #bac, bce, ecf, dbe
290        elements = [ [1,0,2], [1,2,4], [4,2,5], [3,1,4] ]
291
292        domain = Generic_Domain(points, elements)
293        domain.conserved_quantities = ['stage', 'ymomentum']
294        domain.quantities['stage'] =\
295                                   Quantity(domain, [[1,2,3], [5,5,5],
296                                                     [0,0,9], [-6, 3, 3]])
297
298        domain.quantities['ymomentum'] =\
299                                   Quantity(domain, [[2,3,4], [5,5,5],
300                                                     [0,0,9], [-6, 3, 3]])
301
302        domain.check_integrity()
303
304
305        #Write file
306        filename = 'boundarytest' + str(time.time())
307        fid = open(filename + '.txt', 'w')
308        start = time.mktime(time.strptime('2000', '%Y'))
309        dt = 5*60  #Five minute intervals
310        for i in range(10):
311            t = start + i*dt
312            t_string = time.strftime(time_format, time.gmtime(t))
313
314            fid.write('%s,%f %f\n' %(t_string, 1.0*i, sin(i*2*pi/10)))
315        fid.close()
316
317
318        #Convert ASCII file to NetCDF (Which is what we really like!)
319       
320        from anuga.shallow_water.data_manager import timefile2netcdf
321       
322        timefile2netcdf(filename, quantity_names = ['stage', 'ymomentum'])
323       
324
325
326        F = File_boundary(filename + '.tms', domain)
327
328       
329        os.remove(filename + '.txt')
330        os.remove(filename + '.tms')       
331
332
333       
334
335        #Check that midpoint coordinates at boundary are correctly computed
336        assert num.allclose( F.midpoint_coordinates,
337                             [[1.0, 0.0], [0.0, 1.0], [3.0, 0.0],
338                              [3.0, 1.0], [1.0, 3.0], [0.0, 3.0]])
339
340        #assert allclose(F.midpoint_coordinates[(3,2)], [0.0, 3.0])
341        #assert allclose(F.midpoint_coordinates[(3,1)], [1.0, 3.0])
342        #assert allclose(F.midpoint_coordinates[(0,2)], [0.0, 1.0])
343        #assert allclose(F.midpoint_coordinates[(0,0)], [1.0, 0.0])
344        #assert allclose(F.midpoint_coordinates[(2,0)], [3.0, 0.0])
345        #assert allclose(F.midpoint_coordinates[(2,1)], [3.0, 1.0])
346
347
348        #Check time interpolation
349        from anuga.config import default_boundary_tag
350        domain.set_boundary( {default_boundary_tag: F} )
351
352        domain.time = 5*30/2  #A quarter way through first step
353        q = F.evaluate()
354        assert num.allclose(q, [1.0/4, sin(2*pi/10)/4])
355
356
357        domain.time = 2.5*5*60  #Half way between steps 2 and 3
358        q = F.evaluate()
359        assert num.allclose(q, [2.5, (sin(2*2*pi/10) + sin(3*2*pi/10))/2])
360
361
362
363    def test_fileboundary_exception(self):
364        """Test that boundary object complains if number of
365        conserved quantities are wrong
366        """
367
368        from quantity import Quantity
369        import time, os
370        from math import sin, pi
371        from anuga.config import time_format
372
373        a = [0.0, 0.0]
374        b = [0.0, 2.0]
375        c = [2.0,0.0]
376        d = [0.0, 4.0]
377        e = [2.0, 2.0]
378        f = [4.0,0.0]
379
380        points = [a, b, c, d, e, f]
381
382        #bac, bce, ecf, dbe
383        elements = [ [1,0,2], [1,2,4], [4,2,5], [3,1,4] ]
384
385        domain = Generic_Domain(points, elements)
386        domain.conserved_quantities = ['stage', 'xmomentum', 'ymomentum']
387        domain.evolved_quantities = ['stage', 'xmomentum', 'ymomentum']
388        domain.quantities['stage'] =\
389                                   Quantity(domain, [[1,2,3], [5,5,5],
390                                                     [0,0,9], [-6, 3, 3]])
391
392        domain.quantities['xmomentum'] =\
393                                   Quantity(domain, [[2,3,4], [5,5,5],
394                                                     [0,0,9], [-6, 3, 3]])
395        domain.quantities['ymomentum'] =\
396                                   Quantity(domain, [[2,3,4], [5,5,5],
397                                                     [0,0,9], [-6, 3, 3]])
398
399        domain.check_integrity()
400
401        #Write file (with only two values)
402        filename = 'boundarytest' + str(time.time())
403        fid = open(filename + '.txt', 'w')
404        start = time.mktime(time.strptime('2000', '%Y'))
405        dt = 5*60  #Five minute intervals
406        for i in range(10):
407            t = start + i*dt
408            t_string = time.strftime(time_format, time.gmtime(t))
409
410            fid.write('%s,%f %f\n' %(t_string, 1.0*i, sin(i*2*pi/10)))
411        fid.close()
412
413
414        #Convert ASCII file to NetCDF (Which is what we really like!)
415        from anuga.file_conversion.file_conversion import timefile2netcdf
416       
417        timefile2netcdf(filename+'.txt', quantity_names = ['stage', 'xmomentum'])
418
419       
420        try:
421            F = File_boundary(filename + '.tms',
422                              domain)
423        except:
424            pass
425        else:
426            raise Exception('Should have raised an exception')
427       
428        os.remove(filename + '.txt')
429        os.remove(filename + '.tms')       
430
431
432#-------------------------------------------------------------
433
434if __name__ == "__main__":
435    suite = unittest.makeSuite(Test_Generic_Boundary_Conditions, 'test')
436    runner = unittest.TextTestRunner()
437    runner.run(suite)
Note: See TracBrowser for help on using the repository browser.