source: branches/numpy/anuga/shallow_water/test_tsunami_okada.py @ 6304

Last change on this file since 6304 was 6304, checked in by rwilson, 16 years ago

Initial commit of numpy changes. Still a long way to go.

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