# source:anuga_work/development/demos/beach_parallel.py@4539

Last change on this file since 4539 was 4539, checked in by ole, 16 years ago

Moved files from examples to anuga_work

File size: 4.1 KB
Line
1"""Example of the inundationmodel.
2
3A wave of water is washed up ontp a hypothetical beach.
4This one uses the parallel api
5
6To run:
7
8mpirun -c 4 python beach.py
9"""
10
11######################
12# Module imports
13
14from math import pi
15import time
16
17from anuga.shallow_water import Domain
18from anuga.shallow_water import Reflective_boundary
19from anuga.shallow_water import Dirichlet_boundary
20from anuga.shallow_water import Time_boundary
21from anuga.utilities.polygon import Polygon_function
22
23from Numeric import choose, greater, ones, sin, exp
24
25from anuga_parallel.parallel_api import myid, numprocs, distribute
26
27
28
29#------------------
30# Define geometries
31#------------------
32
33def bathymetry(x,y):
34    cut = 75
35    res = choose( greater(x, cut), ( -(x - 55)/5, -4*ones(x.shape) ))
36    res == (100-y)/50 + 1
37    return res
38
39
40def topography(x, y):
41    z = -4.0*x/25 + 8  + (100-y)/50                         # Beach
42    z += 6*exp( -((x-30)/10)**2 ) * exp( -((y-50)/8)**2 )   # Mound
43    z += 4*exp( -((x-30)/4)**2 ) * exp( -((y-26)/10)**2 )   # Extra ridge
44    z += 4*exp( -((x-30)/5)**2 ) * exp( -((y-10)/8)**2 )    # Extra ridge
45    z -= 4*exp( -((x-15)/6)**2 ) * exp( -((y-20)/12)**2 )   # Depression
46    z += 1.2*exp( -((x-88)/7)**2 ) + exp( -((y-20)/25)**2 ) # Seafloor
47    return z
48
49
50def riverbed(x, y):
51    return (y-100)/70 - 4.0*x/25 + 8
52
53
54# Polygons
55shoreline = [[40,0], [100,0], [100,100], [65,100],
56             [55,90], [55,70], [56,50], [50,25], [40,0]]
57
58land = [[65,100], [55,90], [55,70], [56,50], [50,25], [40,0],
59        [0,0], [0,100]]
60
61water = [[55,0], [100,0], [100,100], [55,100]]
62all = [[0,0], [0,100], [100,100], [100,0]]
63
64
65building1 = [[45,80], [49,78], [52,83], [46,83]]
66building2 = [[35,75], [40,75], [40,80], [35,80]]
67building3 = [[42,63], [46,61], [48,65], [44,67]]
68building4 = [[28,56], [28,51], [33,51], [33,56]]
69building5 = [[10,70], [10,65], [15,65], [15,70]]
70building6 = [[10,50], [10,45], [15,45], [15,50]]
71
72river = [[20,100], [18,90], [20,80], [20,60], [15,40], [11,20], [2,0], [10,0],
73         [14,10], [20,30], [24,45], [27,80], [27,85], [35,90], [39,100]]
74
75
76
77#----------------------
78# Domain
79#----------------------
80name = 'beach'
81print 'Creating domain from %s.tsh' %name
82domain = Domain(mesh_filename = name + '.tsh',
83                use_cache=True, verbose=True)
84domain.set_name(name)
85domain.set_default_order(2)
86
87
88#----------------------
89# Initial conditions
90#----------------------
91domain.set_quantity('elevation',
92                    Polygon_function( [(all, topography),
93                                       (building1, 7), (building2, 8),
94                                       (building3, 7), (building4, 13),
95                                       (building5, 10), (building6, 11)]))
96
97domain.set_quantity('stage',
98                    Polygon_function( [(water, -1.5),
99                                       (land, -10)] ))
100
101domain.set_quantity('friction', 0.03)
102print domain.statistics()
103
104
105#----------------------
106# Boundary conditions
107#----------------------
108Br = Reflective_boundary(domain)
109Bd = Dirichlet_boundary([-12, 0.0, 0.0])
110#Bt = Time_boundary(domain, lambda t: [ 3.0*(1+sin(2*pi*t/100)), 0.0, 0.0])
111
112tags = {}
113tags['wall'] = Br
114tags['wall1'] = Br
115tags['outflow'] = Bd
116tags['exterior'] = Br
117tags['external'] = Br
118tags['land'] = Bd
119tags['westbank'] = None    #Riverbank
120tags['eastbank'] = None    #Riverbank
121tags['eastbankN'] = None   #Riverbank
122tags['ocean'] = None       # Bind this one later
123
124domain.set_boundary(tags)
125
126
127#--------------------
128# Distribute domain
129#--------------------
130domain = distribute(domain)
131
132Bt = Time_boundary(domain, lambda t: [ 4.0*(1+sin(2*pi*t/50)), -1.0, 0.0])
133domain.set_boundary({'ocean': Bt})
134
135#----------------------
136# Evolve through time
137#----------------------
138t0 = time.time()
139for t in domain.evolve(yieldstep = 0.2, finaltime = 300):
140    #domain.write_time()
141
142    w = domain.get_maximum_inundation_elevation()
143    x, y = domain.get_maximum_inundation_location()
144    t = domain.get_time()
145    print '  Coastline elevation at t = %.2f is %.2f at loc (x,y)=(%.2f, %.2f)' %(t, w, x, y)
146
147
148print 'Simulation took %.2f seconds' %(time.time()-t0)
149
Note: See TracBrowser for help on using the repository browser.