source: trunk/anuga_core/source/anuga/shallow_water/test_tsunami_okada.py @ 7780

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

Almost all failing tests fixed.

File size: 9.0 KB
Line 
1import unittest
2import numpy as num
3from tsunami_okada import earthquake_tsunami,Okada_func
4from anuga.shallow_water.shallow_water_domain import Domain
5
6class Test_eq(unittest.TestCase):
7    def setUp(self):
8        pass
9
10    def tearDown(self):
11        pass
12
13
14    def test_Okada_func(self):
15        from os import sep, getenv
16        import sys
17        from anuga.abstract_2d_finite_volumes.mesh_factory \
18        import rectangular_cross
19
20        from anuga.abstract_2d_finite_volumes.quantity import Quantity
21        from anuga.utilities.system_tools import get_pathname_from_package
22        """
23        Pick the test you want to do; T= 0 test a point source,
24        T= 1  test single rectangular source, T= 2 test multiple
25        rectangular sources
26        """
27        # Get path where this test is run
28        path = get_pathname_from_package('anuga.shallow_water')
29       
30        # Choose what test to proceed
31        T = 1
32       
33
34        if T==0:
35            # Fortran output file           
36            filename = path+sep+'fullokada_SP.txt'
37           
38            # Initial condition of earthquake for multiple source
39            x0 = 7000.0
40            y0 = 10000.0
41            length = 0
42            width =0
43            strike = 0.0
44            depth = 15.0
45            slip = 10.0
46            dip =15.0
47            rake =90.0
48            ns=1
49            NSMAX=1
50        elif T==1:
51            # Fortran output file       
52            filename = path+sep+'fullokada_SS.txt'
53           
54            # Initial condition of earthquake for multiple source
55            x0 = 7000.0
56            y0 = 10000.0
57            length = 10.0
58            width =6.0
59            strike = 0.0
60            depth = 15.0
61            slip = 10.0
62            dip =15.0
63            rake =90.0
64            ns=1
65            NSMAX=1
66           
67        elif T==2:
68
69            # Fortran output file
70            filename = path+sep+'fullokada_MS.txt'
71           
72            # Initial condition of earthquake for multiple source
73            x0 = [7000.0,10000.0]
74            y0 = [10000.0,7000.0]
75            length = [10.0,10.0]
76            width =[6.0,6.0]
77            strike = [0.0,0.0]
78            depth = [15.0,15.0]
79            slip = [10.0,10.0]
80            dip = [15.0,15.0]
81            rake = [90.0,90.0]
82            ns=2
83            NSMAX=2
84
85
86
87        # Get output file from original okada fortran script.
88        # Vertical displacement is listed under tmp.
89        polyline_file=open(filename,'r')
90        lines=polyline_file.readlines()
91        polyline_file.close()
92        tmp=[]
93        stage=[]
94        for line in lines [0:]:
95            field = line.split('    ')
96            z=float(field[2])
97            tmp.append(z)
98
99       
100
101               
102        #create domain
103        dx = dy = 4000
104        l=100000
105        w=100000
106        #create topography
107        def topography(x,y):
108            el=-1000
109            return el
110       
111        points, vertices, boundary = rectangular_cross(int(l/dx), int(w/dy),
112                                               len1=l, len2=w)
113        domain = Domain(points, vertices, boundary)   
114        domain.set_name('test')
115        domain.set_quantity('elevation',topography)
116       
117        #create variable with elevation data to implement in okada
118        zrec0 = Quantity(domain)
119        zrec0.set_values(0.0)
120        zrec=zrec0.get_vertex_values(xy=True)
121        # call okada
122        Ts= Okada_func(ns=ns, NSMAX=NSMAX,length=length, width=width, dip=dip, \
123                       x0=x0, y0=y0, strike=strike, depth=depth, \
124                       slip=slip, rake=rake,zrec=zrec)
125
126        #create a variable to store vertical displacement throughout the domain
127        tsunami = Quantity(domain)
128        tsunami.set_values(Ts)
129
130        # get vertical displacement at each point of the domain respecting
131        # original script's order
132        interpolation_points=[]
133        k=0.0
134        for i in range(0,6):
135            for j in range(0,6):
136                p=j*4000
137                Yt=p
138                Xt=k
139                interpolation_points.append([Xt, Yt])
140
141            k=k+4000
142        Z=tsunami.get_values(interpolation_points=interpolation_points,
143                             location='edges')
144
145        stage = -Z # FIXME(Ole): Why the sign flip?
146                   # Displacement in fortran code is looking downward
147        #print tmp
148        #print 'hello',stage
149        assert num.allclose(stage,tmp,atol=1.e-3)
150       
151    def test_earthquake_tsunami(self):
152        from os import sep, getenv
153        import sys
154        from anuga.abstract_2d_finite_volumes.mesh_factory \
155                        import rectangular_cross
156        from anuga.abstract_2d_finite_volumes.quantity import Quantity
157        from anuga.utilities.system_tools import get_pathname_from_package
158        """
159        Pick the test you want to do; T= 0 test a point source,
160        T= 1  test single rectangular source, T= 2 test multiple
161        rectangular sources
162        """
163
164        # Get path where this test is run
165        path= get_pathname_from_package('anuga.shallow_water')
166       
167        # Choose what test to proceed
168        T=1
169
170        if T==0:
171            # Fortran output file
172            filename = path+sep+'fullokada_SP.txt'
173           
174            # Initial condition of earthquake for multiple source
175            x0 = 7000.0
176            y0 = 10000.0
177            length = 0
178            width =0
179            strike = 0.0
180            depth = 15.0
181            slip = 10.0
182            dip =15.0
183            rake =90.0
184            ns=1
185            NSMAX=1
186        elif T==1:
187            # Fortran output file
188            filename = path+sep+'fullokada_SS.txt'
189           
190            # Initial condition of earthquake for multiple source
191            x0 = 7000.0
192            y0 = 10000.0
193            length = 10.0
194            width =6.0
195            strike = 0.0
196            depth = 15.0
197            slip = 10.0
198            dip =15.0
199            rake =90.0
200            ns=1
201            NSMAX=1
202           
203        elif T==2:
204
205            # Fortran output file
206            filename = path+sep+'fullokada_MS.txt'
207           
208            # Initial condition of earthquake for multiple source
209            x0 = [7000.0,10000.0]
210            y0 = [10000.0,7000.0]
211            length = [10.0,10.0]
212            width =[6.0,6.0]
213            strike = [0.0,0.0]
214            depth = [15.0,15.0]
215            slip = [10.0,10.0]
216            dip = [15.0,15.0]
217            rake = [90.0,90.0]
218            ns=2
219            NSMAX=2
220
221
222
223        # Get output file from original okada fortran script.
224        # Vertical displacement is listed under tmp.
225        polyline_file=open(filename,'r')
226        lines=polyline_file.readlines()
227        polyline_file.close()
228        tmp=[]
229        stage=[]
230        for line in lines [0:]:
231            field = line.split('    ')
232            z=float(field[2])
233            tmp.append(z)
234
235         
236        # Create domain
237        dx = dy = 4000
238        l=20000
239        w=20000
240       
241        # Create topography
242        def topography(x,y):
243            el=-1000
244            return el
245       
246        points, vertices, boundary = rectangular_cross(int(l/dx), int(w/dy),
247                                               len1=l, len2=w)
248        domain = Domain(points, vertices, boundary)   
249        domain.set_name('test')
250        domain.set_quantity('elevation',topography)
251        Ts = earthquake_tsunami(ns=ns,NSMAX=NSMAX,length=length, width=width, strike=strike,\
252                                depth=depth,dip=dip, xi=x0, yi=y0,z0=0, slip=slip, rake=rake,\
253                                domain=domain, verbose=False)
254       
255        # Create a variable to store vertical displacement throughout the domain
256        tsunami = Quantity(domain)
257        tsunami.set_values(Ts)
258        interpolation_points=[]
259
260        #k=0.0
261        #for i in range(0,6):
262        #    for j in range(0,6):
263        #        p=j*4000
264        #        Yt=p
265        #        Xt=k
266        #        Z=tsunami.get_values(interpolation_points=[[Xt,Yt]]
267        #                             ,location='edges')
268        #        stage.append(-Z[0])
269        #    k=k+4000
270        #
271        #assert allclose(stage,tmp,atol=1.e-3)
272
273        # Here's a faster way - try that in the first test
274        interpolation_points=[]
275        k=0.0
276        for i in range(0,6):
277            for j in range(0,6):
278                p=j*4000
279                Yt=p
280                Xt=k
281                interpolation_points.append([Xt, Yt])
282
283            k=k+4000
284        Z=tsunami.get_values(interpolation_points=interpolation_points,
285                             location='edges')
286
287        stage = -Z # FIXME(Ole): Why the sign flip?
288                   # Displacement in fortran code is looking downward
289        #print 'c est fini'
290        #print tmp
291        #print 'hello',stage   
292        assert num.allclose(stage,tmp,atol=1.e-3)
293
294#-------------------------------------------------------------
295
296if __name__ == "__main__":
297    suite = unittest.makeSuite(Test_eq,'test')
298    runner = unittest.TextTestRunner()
299    runner.run(suite)
300
Note: See TracBrowser for help on using the repository browser.