source: anuga_core/source/anuga/shallow_water/eq_test.py @ 5474

Last change on this file since 5474 was 4436, checked in by nick, 17 years ago

clean up and update source function 'eqf however there are still problems with it

File size: 3.2 KB
Line 
1
2from anuga.shallow_water import Domain, Reflective_boundary, Dirichlet_boundary
3from eqf_v2 import earthquake_tsunami
4
5res = 10000000.
6intres = 1000000.
7boundingpoly=[[-100000,-80000],[-100000,80000],[100000,80000],[100000,-80000]]
8poly = [[-40000,-20000],[40000,-20000],[40000,30000],[-40000,30000]]
9interior_regions = [[poly,intres]]
10
11from anuga.pmesh.mesh_interface import create_mesh_from_regions
12from caching import cache
13meshname = 'another_test'+'.msh'
14
15_ = cache(create_mesh_from_regions,
16          boundingpoly,
17           {'boundary_tags': {'e0': [0], 'e1': [1], 'e2': [2], 'e3': [3]},
18           'maximum_triangle_area': res,
19           'filename': meshname,
20           'interior_regions': interior_regions},
21           verbose = True)
22
23domain = Domain(meshname, use_cache = True, verbose = True)
24
25basename = 'test'
26domain.set_name(basename)
27domain.set_quantities_to_be_stored(['stage', 'xmomentum', 'ymomentum'])
28domain.set_minimum_storable_height(0.01)
29
30tsunami_source = earthquake_tsunami(length=25.0,
31                                    width=5.0,
32                                    strike=0.,
33                                    depth=3.50,
34#                                    depth=3500.0,
35                                    dip=15.0,
36                                    slip=1.,
37                                    x0=0.0,
38                                    y0=0.0,
39                                    domain=domain,
40                                    verbose=True)
41
42stage0 = 1.0
43
44#sets tsunami_source to the domain
45domain.set_quantity('stage', tsunami_source)
46z2 = domain.get_quantity('stage')
47int2 = z2.get_integral()/400000.0
48#print 'integral before and after', int1, int2
49
50print 'len(z2)',len(z2)
51print 'dir(z2)',dir(z2)
52
53domain.set_quantity('elevation', -10.0, alpha = 0.1)
54
55Br = Reflective_boundary(domain)
56Bd = Dirichlet_boundary([0,0,0])
57domain.set_boundary( {'e0': Br,  'e1': Br, 'e2': Br, 'e3': Br} )
58
59import time
60
61t0 = time.time()
62
63for t in domain.evolve(yieldstep = 1, finaltime = 1): 
64    domain.write_time()
65    domain.write_boundary_statistics(tags = 'e2')
66   
67from anuga.abstract_2d_finite_volumes.util import file_function
68from Numeric import array, Float
69from pylab import plot, ion, hold,savefig
70
71'''
72#points=array([[0,0],[0,10]])
73max = 100
74skip = 1000
75
76y = 10000
77points=array([[0]*2]*max)
78print points
79half_max_skip=(max*skip)/2
80for i in range(max):
81    print i*skip-half_max_skip
82    points[i]=[i*skip-half_max_skip, y]
83#    print i
84'''
85interval=500
86profile_lenght= 20000
87number_points = profile_lenght/interval
88y = 10000
89points=array([[0]*2]*number_points)
90print points
91half_profile=profile_lenght/2
92for i in range(number_points):
93    print i*interval-half_profile
94    points[i]=[i*interval-half_profile, y]
95
96
97print "points[0,:]",points
98ion()
99print 'hello'
100
101F = file_function('test.sww', quantities = 'stage', interpolation_points=points,verbose = True,use_cache=True)
102
103print F.statistics()
104t=1
105x = 0.0
106y = 0.0
107#print F(t=1, x=5000, y=5000)
108profile=[]
109
110for i in range(number_points):
111    profile.append(F(0,point_id=i))
112    print i, F(0,point_id=i)
113print'profile', profile#,points[:,1]
114
115plot(points[:,0],profile)
116savefig("profile",dpi=300)
117
118   
119   
120
121
Note: See TracBrowser for help on using the repository browser.