1 | """Script for running a dam break simulation of UQ's dam break tank. |
---|
2 | |
---|
3 | |
---|
4 | Ole Nielsen and Duncan Gray, GA - 2006 |
---|
5 | """ |
---|
6 | |
---|
7 | |
---|
8 | #---------------------------------------------------------------------------- |
---|
9 | # Import necessary modules |
---|
10 | #---------------------------------------------------------------------------- |
---|
11 | |
---|
12 | # Standard modules |
---|
13 | import time |
---|
14 | import sys |
---|
15 | from shutil import copy |
---|
16 | from os import path |
---|
17 | |
---|
18 | # Related major packages |
---|
19 | from anuga.shallow_water import Domain, Reflective_boundary, \ |
---|
20 | Dirichlet_boundary, Time_boundary, File_boundary |
---|
21 | from anuga.abstract_2d_finite_volumes.util import start_screen_catcher, \ |
---|
22 | copy_code_files |
---|
23 | from anuga.abstract_2d_finite_volumes.region import Set_region |
---|
24 | from anuga.fit_interpolate.interpolate import interpolate_sww2csv |
---|
25 | import create_mesh |
---|
26 | #from anuga.pyvolution.data_manager import timefile2netcdf |
---|
27 | from math import cos |
---|
28 | from math import cosh |
---|
29 | from math import sin |
---|
30 | from math import sqrt |
---|
31 | from math import fabs |
---|
32 | |
---|
33 | |
---|
34 | |
---|
35 | # Application specific imports |
---|
36 | import project # Definition of file names and polygons |
---|
37 | |
---|
38 | |
---|
39 | def main(): |
---|
40 | |
---|
41 | #solitary wave solution parameters |
---|
42 | slopes = [0.1] #plane beach slope (0 if bathymetry is read from file) |
---|
43 | |
---|
44 | H_d_ratios = [0.00092] #ratio wave height to initial depth |
---|
45 | |
---|
46 | # pos = [0.002, 0.02, 0.08, 0.32] |
---|
47 | |
---|
48 | for slope in slopes: |
---|
49 | for H_d_ratio in H_d_ratios: |
---|
50 | # for po in pos: |
---|
51 | scenario(slope, H_d_ratio) |
---|
52 | |
---|
53 | def scenario(slope, H_d_ratio): |
---|
54 | |
---|
55 | longshore_length = 150 |
---|
56 | offshore_depth = 98 |
---|
57 | offshore_length = offshore_depth/slope |
---|
58 | bay_length = offshore_depth/4/slope |
---|
59 | period = 20 |
---|
60 | x1 = 4000 |
---|
61 | |
---|
62 | #------------------------------------------------------------------------- |
---|
63 | # Setup archiving of simulations |
---|
64 | #------------------------------------------------------------------------- |
---|
65 | |
---|
66 | |
---|
67 | id = 'Sine_wave_plane_beach'+'_s'+ str(slope)+ '_wr'+ str(H_d_ratio) |
---|
68 | copy (project.codedirname, project.outputtimedir + 'project.py') |
---|
69 | run_name = 'run_dam.py' |
---|
70 | run_name_out = 'run_dam'+id+'.py' |
---|
71 | copy (project.codedir + 'run_dam.py', project.outputtimedir + run_name_out) |
---|
72 | copy (project.codedir + 'create_mesh.py', |
---|
73 | project.outputtimedir + 'create_mesh.py') |
---|
74 | print'output dir', project.outputtimedir |
---|
75 | |
---|
76 | #FIXME this isn't working |
---|
77 | #normal screen output is stored in |
---|
78 | screen_output_name = project.outputtimedir + "screen_output.txt" |
---|
79 | screen_error_name = project.outputtimedir + "screen_error.txt" |
---|
80 | |
---|
81 | #------------------------------------------------------------------------- |
---|
82 | # Create the triangular mesh |
---|
83 | #------------------------------------------------------------------------- |
---|
84 | |
---|
85 | create_mesh.generate(project.mesh_filename, |
---|
86 | offshore_length, bay_length, longshore_length, |
---|
87 | is_course=True) # this creates the mesh |
---|
88 | #is_course=False) # this creates the mesh |
---|
89 | |
---|
90 | head,tail = path.split(project.mesh_filename) |
---|
91 | copy (project.mesh_filename, |
---|
92 | project.outputtimedir + tail ) |
---|
93 | #------------------------------------------------------------------------- |
---|
94 | # Setup computational domain |
---|
95 | #------------------------------------------------------------------------- |
---|
96 | domain = Domain(project.mesh_filename, use_cache = False, verbose = True) |
---|
97 | |
---|
98 | |
---|
99 | print 'Number of triangles = ', len(domain) |
---|
100 | |
---|
101 | print 'The extent is ', domain.get_extent() |
---|
102 | print domain.statistics() |
---|
103 | |
---|
104 | |
---|
105 | domain.set_name(project.basename + id) |
---|
106 | domain.set_datadir(project.outputtimedir) |
---|
107 | domain.set_quantities_to_be_stored(['stage', 'xmomentum', 'ymomentum']) |
---|
108 | domain.set_minimum_storable_height(0.01) |
---|
109 | domain.set_store_vertices_uniquely(True) # for writing to sww |
---|
110 | |
---|
111 | #------------------------------------------------------------------------- |
---|
112 | # Setup initial conditions |
---|
113 | #------------------------------------------------------------------------- |
---|
114 | g = 9.81 |
---|
115 | |
---|
116 | def elevation_tilt(x, y): |
---|
117 | return x*slope |
---|
118 | |
---|
119 | domain.set_quantity('stage', filename = 'stage.xya', |
---|
120 | alpha = 0.001,verbose = True, use_cache = True) |
---|
121 | domain.set_quantity('friction', 0) |
---|
122 | domain.set_quantity('elevation', elevation_tilt) |
---|
123 | |
---|
124 | #bathymetry read from file |
---|
125 | #domain.set_quantity('elevation', filename = 'elevation.xya', |
---|
126 | # alpha = 10.0,verbose = True, use_cache = True) |
---|
127 | |
---|
128 | print 'Available boundary tags', domain.get_boundary_tags() |
---|
129 | |
---|
130 | Br = Reflective_boundary(domain) |
---|
131 | |
---|
132 | #bore described by function (solitary wave solution) |
---|
133 | W = Time_boundary(domain = domain, |
---|
134 | |
---|
135 | f=lambda t: [-offshore_depth*H_d_ratio*cos(2*3.14*t/period), 0, 0]) #sine wave |
---|
136 | |
---|
137 | # f=lambda t: [offshore_depth*H_d_ratio/cosh(sqrt(3*H_d_ratio*po/4)* |
---|
138 | # (sqrt(g/offshore_depth)*t-x1/offshore_depth))/cosh(sqrt(3*H_d_ratio*po/4)* |
---|
139 | # (sqrt(g/offshore_depth)*t-x1/offshore_depth)), 0, 0]) |
---|
140 | |
---|
141 | domain.set_boundary( {'wall': Br, 'edge': Br, 'back':W} ) |
---|
142 | |
---|
143 | #------------------------------------------------------------------------- |
---|
144 | # Evolve system through time |
---|
145 | #------------------------------------------------------------------------- |
---|
146 | t0 = time.time() |
---|
147 | |
---|
148 | for t in domain.evolve(yieldstep = 1, finaltime =200): #enter timestep and final time |
---|
149 | domain.write_time() |
---|
150 | |
---|
151 | print 'That took %.2f seconds' %(time.time()-t0) |
---|
152 | print 'finished' |
---|
153 | |
---|
154 | points = [[-offshore_length, longshore_length / 2], |
---|
155 | [(-offshore_length)*0.95, longshore_length / 2], |
---|
156 | [(-offshore_length)*0.9, longshore_length / 2], |
---|
157 | [(-offshore_length)*0.85, longshore_length / 2], |
---|
158 | [(-offshore_length)*0.8, longshore_length / 2], |
---|
159 | [(-offshore_length)*0.75, longshore_length / 2], |
---|
160 | [(-offshore_length)*0.7, longshore_length / 2], |
---|
161 | [(-offshore_length)*0.65, longshore_length / 2], |
---|
162 | [(-offshore_length)*0.6, longshore_length / 2], |
---|
163 | [(-offshore_length)*0.55, longshore_length / 2], |
---|
164 | [(-offshore_length)*0.5, longshore_length / 2], |
---|
165 | [(-offshore_length)*0.45, longshore_length / 2], |
---|
166 | [(-offshore_length)*0.4, longshore_length / 2], |
---|
167 | [(-offshore_length)*0.35, longshore_length / 2], |
---|
168 | [(-offshore_length)*0.3, longshore_length / 2], |
---|
169 | [(-offshore_length)*0.25, longshore_length / 2], |
---|
170 | [(-offshore_length)*0.2, longshore_length / 2], |
---|
171 | [(-offshore_length)*0.15, longshore_length / 2], |
---|
172 | [(-offshore_length)*0.1, longshore_length / 2], |
---|
173 | [(-offshore_length)*0.09, longshore_length / 2], |
---|
174 | [(-offshore_length)*0.08, longshore_length / 2], |
---|
175 | [(-offshore_length)*0.07, longshore_length / 2], |
---|
176 | [(-offshore_length)*0.06, longshore_length / 2], |
---|
177 | [(-offshore_length)*0.05, longshore_length / 2], |
---|
178 | [(-offshore_length)*0.04, longshore_length / 2], |
---|
179 | [(-offshore_length)*0.03, longshore_length / 2], |
---|
180 | [(-offshore_length)*0.02, longshore_length / 2], |
---|
181 | [(-offshore_length)*0.01, longshore_length / 2], |
---|
182 | [(-offshore_length)*0, longshore_length / 2], |
---|
183 | [(offshore_length)*0.01, longshore_length / 2], |
---|
184 | [(offshore_length)*0.02, longshore_length / 2], |
---|
185 | [(offshore_length)*0.02, longshore_length / 2], |
---|
186 | [(offshore_length)*0.03, longshore_length / 2], |
---|
187 | [(offshore_length)*0.04, longshore_length / 2]] |
---|
188 | #------------------------------------------------------------------------- |
---|
189 | # Calculate gauge info |
---|
190 | #------------------------------------------------------------------------- |
---|
191 | interpolate_sww2csv(project.outputtimedir + project.basename \ |
---|
192 | + id+".sww", |
---|
193 | points, |
---|
194 | project.depth_filename + id + '.csv', |
---|
195 | project.velocity_x_filename + id + '.csv', |
---|
196 | project.velocity_y_filename + id + '.csv') |
---|
197 | |
---|
198 | #------------------------------------------------------------- |
---|
199 | if __name__ == "__main__": |
---|
200 | main() |
---|
201 | |
---|