source: anuga_validation/uq_flume_2007/run_dam.py @ 4773

Last change on this file since 4773 was 4773, checked in by duncan, 15 years ago

Working on the anuga-2007 paper. Adding results from the UQ flume.

File size: 4.1 KB
Line 
1"""Script for running a dam break simulation of UQ's dam break tank.
2
3The simulation is based on the simulation Matt has presented and will
4form part of the ANUGA validation paper.
5
6i want to compare the stage at .4m and the velocity at .45m.
7
8
9Ole Nielsen and Duncan Gray, GA - 2006
10"""
11
12
13#----------------------------------------------------------------------------
14# Import necessary modules
15#----------------------------------------------------------------------------
16
17# Standard modules
18import time
19import sys
20from shutil import copy
21from os import path
22
23# Related major packages
24from anuga.shallow_water import Domain, Reflective_boundary, \
25                            Dirichlet_boundary, Time_boundary, File_boundary
26from anuga.abstract_2d_finite_volumes.region import Set_region
27from anuga.fit_interpolate.interpolate import interpolate_sww2csv
28from anuga.abstract_2d_finite_volumes.util import start_screen_catcher, \
29     copy_code_files
30import create_mesh
31
32import project
33
34
35# Application specific imports
36
37
38def main():
39     
40    slope= 0
41    friction = 0.01 
42    inital_depth = 0.2
43    gate_position = 0.75
44   
45    return scenario(slope, friction, inital_depth, gate_position)
46           
47
48def scenario(slope, friction, inital_depth, gate_position):
49
50    #-------------------------------------------------------------------------
51    # Create the triangular mesh
52    #-------------------------------------------------------------------------
53
54    create_mesh.generate(project.mesh_filename,
55                         gate_position,
56                         #is_course=True) # this creates the mesh
57                         is_course=False) # this creates the mesh
58
59    head,tail = path.split(project.mesh_filename)
60    #-------------------------------------------------------------------------
61    # Setup computational domain
62    #-------------------------------------------------------------------------
63    domain = Domain(project.mesh_filename, use_cache = False, verbose = True)
64   
65
66    print 'Number of triangles = ', len(domain)
67    print 'The extent is ', domain.get_extent()
68    print domain.statistics()
69
70   
71    domain.set_name(project.basename)
72    domain.set_datadir('.')
73    domain.set_quantities_to_be_stored(['stage', 'xmomentum', 'ymomentum'])
74    domain.set_minimum_storable_height(0.001)
75    domain.set_store_vertices_uniquely(True)  # for writting to sww
76
77    #-------------------------------------------------------------------------
78    # Setup initial conditions
79    #-------------------------------------------------------------------------
80
81
82    def elevation_tilt(x, y):
83        return x*slope
84       
85    domain.set_quantity('stage', elevation_tilt)
86    domain.set_quantity('friction', friction) 
87    domain.set_quantity('elevation',elevation_tilt)
88
89    print 'Available boundary tags', domain.get_boundary_tags()
90    domain.set_region('dam','stage',inital_depth,
91                                 location = 'unique vertices') 
92
93    Br = Reflective_boundary(domain)
94    Bd = Dirichlet_boundary([0,0,0])  # to drain the water out.
95    domain.set_boundary( {'wall': Br, 'edge': Bd} )
96
97    #-------------------------------------------------------------------------
98    # Evolve system through time
99    #-------------------------------------------------------------------------
100    t0 = time.time()
101
102    for t in domain.evolve(yieldstep = 0.1, finaltime = 10): #enter timestep and final time
103        domain.write_time()
104   
105        print 'That took %.2f seconds' %(time.time()-t0)
106        print 'finished'
107       
108    points = [[0.4,0.2], [0.45,0.2]]
109    #-------------------------------------------------------------------------
110    # Calculate gauge info
111    #-------------------------------------------------------------------------
112    print "name",project.basename +".sww"
113    interpolate_sww2csv( project.basename +".sww",
114                        points,
115                        project.depth_filename ,
116                        project.velocity_x_filename,
117                        project.velocity_y_filename )
118
119
120#-------------------------------------------------------------
121if __name__ == "__main__":
122    main()
123   
Note: See TracBrowser for help on using the repository browser.