source: trunk/anuga_work/development/carrier_greenspan/paul_guard/sine_wave.py @ 7924

Last change on this file since 7924 was 5649, checked in by duncan, 17 years ago

A start in the carrier_greenspan validation. Working on it here initially so it doesn't get released until it works.
File size: 8.1 KB
Line 
1"""Script for running a dam break simulation of UQ's dam break tank.
2
3
4Ole Nielsen and Duncan Gray, GA - 2006
5"""
6
7
8#----------------------------------------------------------------------------
9# Import necessary modules
10#----------------------------------------------------------------------------
11
12# Standard modules
13import time
14import sys
15from shutil import copy
16from os import path
17
18# Related major packages
19from anuga.shallow_water import Domain, Reflective_boundary, \
20                            Dirichlet_boundary, Time_boundary, File_boundary
21from anuga.abstract_2d_finite_volumes.util import start_screen_catcher, \
22     copy_code_files
23from anuga.abstract_2d_finite_volumes.region import Set_region
24from anuga.fit_interpolate.interpolate import interpolate_sww2csv
25import create_mesh
26#from anuga.pyvolution.data_manager import timefile2netcdf
27from math import cos
28from math import cosh
29from math import sin
30from math import sqrt
31from math import fabs
32
33
34
35# Application specific imports
36import project                 # Definition of file names and polygons
37
38
39def 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                   
53def 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#-------------------------------------------------------------
199if __name__ == "__main__":
200    main()
201   
Note: See TracBrowser for help on using the repository browser.