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

Last change on this file since 5589 was 5421, checked in by ole, 15 years ago

Changed timestepping_statistics to output real time seconds as well
Also comments and formatting

File size: 8.9 KB
RevLine 
[5204]1import unittest
2from Numeric import allclose
[5244]3from tsunami_okada import earthquake_tsunami,Okada_func
[5204]4
5class Test_eq(unittest.TestCase):
6    def setUp(self):
7        pass
8
9    def tearDown(self):
10        pass
11
12
[5286]13    def test_Okada_func(self):
[5204]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
[5206]20        from anuga.utilities.system_tools import get_pathname_from_package
[5204]21        """
[5206]22        Pick the test you want to do; T= 0 test a point source,
[5204]23        T= 1  test single rectangular source, T= 2 test multiple
24        rectangular sources
25        """
[5206]26        #get path where this test is run
[5421]27        path = get_pathname_from_package('anuga.shallow_water')
[5204]28        #choose what test to proceed
[5421]29        T = 1
[5206]30       
[5204]31
32
33        if T==0:
34            # fotran output file
[5206]35           
36            filename = path+sep+'fullokada_SP.txt'
[5204]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
[5244]43            depth = 15.0
[5204]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
[5206]51            filename = path+sep+'fullokada_SS.txt'
[5204]52        # initial condition of earthquake for multiple source
[5310]53            x0 = 7000.0
54            y0 = 10000.0
55            length = 10.0
56            width =6.0
[5204]57            strike = 0.0
[5310]58            depth = 15.0
59            slip = 10.0
[5204]60            dip =15.0
61            rake =90.0
62            ns=1
63            NSMAX=1
64           
65        elif T==2:
66
67        # fotran output file
[5206]68            filename = path+sep+'fullokada_MS.txt'
[5204]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]
[5244]75            depth = [15.0,15.0]
[5204]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
[5252]101        l=100000
102        w=100000
[5204]103        #create topography
104        def topography(x,y):
[5252]105            el=-1000
[5204]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
[5239]115        zrec0 = Quantity(domain)
[5252]116        zrec0.set_values(0.0)
[5239]117        zrec=zrec0.get_vertex_values(xy=True)
[5204]118        # call okada
[5206]119        Ts= Okada_func(ns=ns, NSMAX=NSMAX,length=length, width=width, dip=dip, \
[5204]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
[5309]129        interpolation_points=[]
[5204]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
[5309]136                interpolation_points.append([Xt, Yt])
137
[5204]138            k=k+4000
[5309]139        Z=tsunami.get_values(interpolation_points=interpolation_points,
140                             location='edges')
[5310]141
142        stage = -Z # FIXME(Ole): Why the sign flip?
143                   # Displacement in fortran code is looking downward
[5312]144        #print tmp
145        #print 'hello',stage
[5310]146        assert allclose(stage,tmp,atol=1.e-3)
[5244]147       
[5204]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
[5206]155        from anuga.utilities.system_tools import get_pathname_from_package
[5204]156        """
[5244]157        Pick the test you want to do; T= 0 test a point source,
[5204]158        T= 1  test single rectangular source, T= 2 test multiple
159        rectangular sources
160        """
[5309]161
[5206]162        #get path where this test is run
163        path= get_pathname_from_package('anuga.shallow_water')
164       
[5204]165        #choose what test to proceed
[5206]166        T=1
[5204]167
168        if T==0:
[5206]169        # fotran output file
170           
171            filename = path+sep+'fullokada_SP.txt'
[5204]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
[5244]178            depth = 15.0
[5204]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
[5206]186            filename = path+sep+'fullokada_SS.txt'
[5204]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
[5244]193            depth = 15.0
[5204]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
[5206]203            filename = path+sep+'fullokada_MS.txt'
[5204]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]
[5244]210            depth = [15.0,15.0]
[5204]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         
[5309]232        # Create domain
[5204]233        dx = dy = 4000
234        l=20000
235        w=20000
[5309]236       
237        # Create topography
[5204]238        def topography(x,y):
[5252]239            el=-1000
[5204]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)
[5206]247        Ts = earthquake_tsunami(ns=ns,NSMAX=NSMAX,length=length, width=width, strike=strike,\
[5252]248                                depth=depth,dip=dip, xi=x0, yi=y0,z0=0, slip=slip, rake=rake,\
[5275]249                                domain=domain, verbose=False)
[5206]250       
[5309]251        # Create a variable to store vertical displacement throughout the domain
[5206]252        tsunami = Quantity(domain)
253        tsunami.set_values(Ts)
[5281]254        interpolation_points=[]
[5309]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=[]
[5204]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
[5281]277                interpolation_points.append([Xt, Yt])
278
[5204]279            k=k+4000
[5281]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
[5312]285        #print 'c est fini'
[5281]286        #print tmp
287        #print 'hello',stage   
288        assert allclose(stage,tmp,atol=1.e-3)
289
[5204]290#-------------------------------------------------------------
291if __name__ == "__main__":
[5309]292    suite = unittest.makeSuite(Test_eq,'test')
[5204]293    runner = unittest.TextTestRunner()
294    runner.run(suite)
295
Note: See TracBrowser for help on using the repository browser.