source: anuga_core/source/anuga/shallow_water/test_tsunami_okada.py @ 5283

Last change on this file since 5283 was 5283, checked in by herve, 17 years ago

changes made according to Ole's suggestion reducing considerably computation time.

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