1 | """Simple water flow example using ANUGA |
---|
2 | |
---|
3 | Water flowing down a channel with more complex topography |
---|
4 | """ |
---|
5 | |
---|
6 | #------------------------------------------------------------------------------ |
---|
7 | # Import necessary modules |
---|
8 | #------------------------------------------------------------------------------ |
---|
9 | from anuga.abstract_2d_finite_volumes.mesh_factory import rectangular_cross |
---|
10 | from anuga.shallow_water import Domain |
---|
11 | from anuga.shallow_water import Reflective_boundary |
---|
12 | from anuga.shallow_water import Dirichlet_boundary |
---|
13 | from anuga.shallow_water import Time_boundary |
---|
14 | from anuga.pmesh.mesh_interface import create_mesh_from_regions |
---|
15 | from anuga.shallow_water import Transmissive_Momentum_Set_Stage_boundary |
---|
16 | from math import tan, sqrt, sin, pi |
---|
17 | |
---|
18 | #------------------------------------------------------------------------------ |
---|
19 | # Project file |
---|
20 | #------------------------------------------------------------------------------ |
---|
21 | #Set up file structure |
---|
22 | anuga_dir = home+'anuga_validation'+sep+'circular_island_tsunami_benchmark'+sep+'anuga'+sep |
---|
23 | meshes_dir = anuga_dir+'meshes'+sep |
---|
24 | meshname = 'circular_mesh.msh' |
---|
25 | output_dir = anuga_dir+'output'+sep |
---|
26 | |
---|
27 | |
---|
28 | ## #------------------------------------------------------------------------------ |
---|
29 | ## # Copy scripts to time stamped output directory and capture screen |
---|
30 | ## # output to file |
---|
31 | ## #------------------------------------------------------------------------------ |
---|
32 | ## print "Processor Name:",get_processor_name() |
---|
33 | ## |
---|
34 | ## #copy script must be before screen_catcher |
---|
35 | ## #print kwargs |
---|
36 | ## |
---|
37 | ## print 'output_dir', output_dir |
---|
38 | ## if myid == 0: |
---|
39 | ## copy_code_files(kwargs['output_dir'],__file__, |
---|
40 | ## dirname(project.__file__)+sep+ project.__name__+'.py' ) |
---|
41 | ## |
---|
42 | ## store_parameters(**kwargs) |
---|
43 | ## |
---|
44 | ## barrier() |
---|
45 | ## |
---|
46 | ## start_screen_catcher(kwargs['output_dir'], myid, numprocs) |
---|
47 | ## |
---|
48 | ## print "Processor Name:",get_processor_name() |
---|
49 | ## |
---|
50 | ## # filenames |
---|
51 | ### meshes_dir_name = project.meshes_dir_name+'.msh' |
---|
52 | ## |
---|
53 | ## # creates copy of code in output dir |
---|
54 | ## print 'min triangles', project.trigs_min, |
---|
55 | ## print 'Note: This is generally about 20% less than the final amount' |
---|
56 | |
---|
57 | |
---|
58 | #------------------------------------------------------------------------------ |
---|
59 | # Setup computational domain |
---|
60 | #------------------------------------------------------------------------------ |
---|
61 | |
---|
62 | |
---|
63 | length = 30. |
---|
64 | width = 25. |
---|
65 | Cx = 12.96 # centre of island on the x axis |
---|
66 | Cy = 13.8 # centre of island on the y axis |
---|
67 | dx = dy = 1 # Resolution: Length of subdivisions on both axes |
---|
68 | water_depth = 0.32 # Can be 0.32 or 0.42 |
---|
69 | |
---|
70 | #boundary |
---|
71 | poly_domain = [[0,0],[length,0],[length,width],[0,width]] |
---|
72 | |
---|
73 | |
---|
74 | # exporting asc grid |
---|
75 | xmin = 0 |
---|
76 | xmax = length |
---|
77 | ymin = 0 |
---|
78 | ymax = width |
---|
79 | |
---|
80 | #Create interior region |
---|
81 | Dome = [[(Cx)-4,(Cy)-4],[(Cx)+4,(Cy)-4,], |
---|
82 | [(Cx)+4,(Cy)+4],[(Cx)-4,(Cy)+4]] |
---|
83 | |
---|
84 | remainder_res=1 |
---|
85 | Dome_res = .01 |
---|
86 | |
---|
87 | interior_dome = [[Dome, Dome_res]] |
---|
88 | |
---|
89 | #Create mesh |
---|
90 | |
---|
91 | print 'start create mesh from regions' |
---|
92 | |
---|
93 | create_mesh_from_regions(poly_domain, |
---|
94 | boundary_tags={'wavemaker': [0], 'right': [1], |
---|
95 | 'top': [2], 'left': [3]}, |
---|
96 | maximum_triangle_area = remainder_res, |
---|
97 | filename=meshname, |
---|
98 | interior_regions = interior_dome, |
---|
99 | use_cache=False, |
---|
100 | verbose=False) |
---|
101 | |
---|
102 | # Setup computational domain |
---|
103 | |
---|
104 | domain = Domain(meshes_dir+sep+meshname, use_cache=False, verbose = True) |
---|
105 | domain.set_name('circular') # Output name |
---|
106 | print 'memory usage before del domain',mem_usage() |
---|
107 | |
---|
108 | print domain.statistics() |
---|
109 | print 'triangles',len(domain) |
---|
110 | |
---|
111 | |
---|
112 | #------------------------------------------------------------------------------ |
---|
113 | # Setup initial conditions |
---|
114 | #------------------------------------------------------------------------------ |
---|
115 | def topography(x,y): |
---|
116 | """Complex topography defined by a function of vectors x and y |
---|
117 | """ |
---|
118 | |
---|
119 | |
---|
120 | z= 0*x # defining z for all values other than the if statements |
---|
121 | r= 3.6 # radius, provided in document |
---|
122 | angle = 14 # angle, provided in document |
---|
123 | h= r*tan(angle/57.2957795) # finding height of cone if not truncated |
---|
124 | |
---|
125 | |
---|
126 | N = len(x) |
---|
127 | for i in range(N): |
---|
128 | |
---|
129 | #truncated top |
---|
130 | if (x[i]-Cx)**2 + (y[i]-Cy)**2 <1.1**2: |
---|
131 | z[i] += 0.625 |
---|
132 | |
---|
133 | # cone |
---|
134 | if (x[i]-Cx)**2 + (y[i]-Cy)**2 <r**2 and (x[i]-Cx)**2 + (y[i]-Cy)**2 >1.1**2: |
---|
135 | z[i] = -(sqrt(((x[i]-Cx)**2+(y[i]-Cy)**2)/((r/h)**2))-h) |
---|
136 | return z |
---|
137 | |
---|
138 | domain.set_quantity('elevation', topography, verbose=True) # Use function for elevation |
---|
139 | domain.set_quantity('friction', 0.01) # Constant friction |
---|
140 | domain.set_quantity('stage',water_depth) # Dry initial condition |
---|
141 | |
---|
142 | |
---|
143 | #------------------------------------------------------------------------------ |
---|
144 | # Setup boundary conditions |
---|
145 | #------------------------------------------------------------------------------ |
---|
146 | # Create boundary function from timeseries provided in file |
---|
147 | |
---|
148 | ##boundary_filename="ts2cnew1_input_20_80sec_new" |
---|
149 | ##prepare_timeboundary(boundary_filename+'.txt') |
---|
150 | ## |
---|
151 | ##function = file_function(boundary_filename+'.tms', |
---|
152 | ## domain, verbose=True) |
---|
153 | def wave_form(t): |
---|
154 | return 0.1*sin(2*pi*t/50.) |
---|
155 | |
---|
156 | # Create and assign boundary objects |
---|
157 | Bw = Dirichlet_boundary([water_depth, 0, 0]) #wall |
---|
158 | Bt = Transmissive_Momentum_Set_Stage_boundary(domain, wave_form) #wavemaker |
---|
159 | |
---|
160 | domain.set_boundary({'left': Bw, 'right': Bw, 'top': Bw, 'wavemaker': Bt}) |
---|
161 | |
---|
162 | |
---|
163 | |
---|
164 | |
---|
165 | #------------------------------------------------------------------------------ |
---|
166 | # Evolve system through time |
---|
167 | #------------------------------------------------------------------------------ |
---|
168 | for t in domain.evolve(yieldstep = 0.2, finaltime = 100.0): |
---|
169 | print domain.timestepping_statistics() |
---|
170 | |
---|
171 | |
---|
172 | |
---|