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

Last change on this file since 5312 was 5312, checked in by herve, 16 years ago

remove print

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 test_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 = 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        # 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        stage = -Z # FIXME(Ole): Why the sign flip?
143                   # Displacement in fortran code is looking downward
144        #print tmp
145        #print 'hello',stage
146        assert allclose(stage,tmp,atol=1.e-3)
147       
148    def test_earthquake_tsunami(self):
149        from os import sep, getenv
150        from Numeric import zeros, Float,allclose
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        # fotran output file
170           
171            filename = path+sep+'fullokada_SP.txt'
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        # fotran output file
186            filename = path+sep+'fullokada_SS.txt'
187        # initial condition of earthquake for multiple source
188            x0 = 7000.0
189            y0 = 10000.0
190            length = 10.0
191            width =6.0
192            strike = 0.0
193            depth = 15.0
194            slip = 10.0
195            dip =15.0
196            rake =90.0
197            ns=1
198            NSMAX=1
199           
200        elif T==2:
201
202        # fotran output file
203            filename = path+sep+'fullokada_MS.txt'
204        # initial condition of earthquake for multiple source
205            x0 = [7000.0,10000.0]
206            y0 = [10000.0,7000.0]
207            length = [10.0,10.0]
208            width =[6.0,6.0]
209            strike = [0.0,0.0]
210            depth = [15.0,15.0]
211            slip = [10.0,10.0]
212            dip = [15.0,15.0]
213            rake = [90.0,90.0]
214            ns=2
215            NSMAX=2
216
217
218
219        # get output file from original okada fortran script.
220        # vertical displacement is listed under tmp.
221        polyline_file=open(filename,'r')
222        lines=polyline_file.readlines()
223        polyline_file.close()
224        tmp=[]
225        stage=[]
226        for line in lines [0:]:
227            field = line.split('    ')
228            z=float(field[2])
229            tmp.append(z)
230
231         
232        # Create domain
233        dx = dy = 4000
234        l=20000
235        w=20000
236       
237        # Create topography
238        def topography(x,y):
239            el=-1000
240            return el
241       
242        points, vertices, boundary = rectangular_cross(int(l/dx), int(w/dy),
243                                               len1=l, len2=w)
244        domain = Domain(points, vertices, boundary)   
245        domain.set_name('test')
246        domain.set_quantity('elevation',topography)
247        Ts = earthquake_tsunami(ns=ns,NSMAX=NSMAX,length=length, width=width, strike=strike,\
248                                depth=depth,dip=dip, xi=x0, yi=y0,z0=0, slip=slip, rake=rake,\
249                                domain=domain, verbose=False)
250       
251        # Create a variable to store vertical displacement throughout the domain
252        tsunami = Quantity(domain)
253        tsunami.set_values(Ts)
254        interpolation_points=[]
255
256        #k=0.0
257        #for i in range(0,6):
258        #    for j in range(0,6):
259        #        p=j*4000
260        #        Yt=p
261        #        Xt=k
262        #        Z=tsunami.get_values(interpolation_points=[[Xt,Yt]]
263        #                             ,location='edges')
264        #        stage.append(-Z[0])
265        #    k=k+4000
266        #
267        #assert allclose(stage,tmp,atol=1.e-3)
268
269        # Here's a faster way - try that in the first test
270        interpolation_points=[]
271        k=0.0
272        for i in range(0,6):
273            for j in range(0,6):
274                p=j*4000
275                Yt=p
276                Xt=k
277                interpolation_points.append([Xt, Yt])
278
279            k=k+4000
280        Z=tsunami.get_values(interpolation_points=interpolation_points,
281                             location='edges')
282
283        stage = -Z # FIXME(Ole): Why the sign flip?
284                   # Displacement in fortran code is looking downward
285        #print 'c est fini'
286        #print tmp
287        #print 'hello',stage   
288        assert allclose(stage,tmp,atol=1.e-3)
289
290#-------------------------------------------------------------
291if __name__ == "__main__":
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.