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

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

Added sww module.

File size: 5.6 KB
Line 
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
10from sww import load_sww_as_domain
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
160#################################################################################
161
162if __name__ == "__main__":
163    suite = unittest.makeSuite(Test_sww, 'test')
164    runner = unittest.TextTestRunner(verbosity=1)
165    runner.run(suite)
Note: See TracBrowser for help on using the repository browser.