source: anuga_validation/uq_flume_2007/run_dam.py @ 7773

Last change on this file since 7773 was 7773, checked in by hudson, 13 years ago

Updated flume validation to use new API.

File size: 4.0 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.utilities.file_utils import copy_code_files
29import create_mesh
30
31import project
32
33
34# Application specific imports
35
36
37def main():
38     
39    slope= 0
40    friction = 0.01 
41    inital_depth = 0.2
42    gate_position = 0.75
43   
44    return scenario(slope, friction, inital_depth, gate_position)
45           
46
47def scenario(slope, friction, inital_depth, gate_position):
48
49    #-------------------------------------------------------------------------
50    # Create the triangular mesh
51    #-------------------------------------------------------------------------
52
53    create_mesh.generate(project.mesh_filename,
54                         gate_position,
55                         #is_course=True) # this creates the mesh
56                         is_course=False) # this creates the mesh
57
58    head,tail = path.split(project.mesh_filename)
59    #-------------------------------------------------------------------------
60    # Setup computational domain
61    #-------------------------------------------------------------------------
62    domain = Domain(project.mesh_filename, use_cache = False, verbose = True)
63   
64
65    print 'Number of triangles = ', len(domain)
66    print 'The extent is ', domain.get_extent()
67    print domain.statistics()
68
69   
70    domain.set_name(project.basename)
71    domain.set_datadir('.')
72    domain.set_quantities_to_be_stored(['stage', 'xmomentum', 'ymomentum'])
73    domain.set_minimum_storable_height(0.001)
74    domain.set_store_vertices_uniquely(True)  # for writting to sww
75
76    #-------------------------------------------------------------------------
77    # Setup initial conditions
78    #-------------------------------------------------------------------------
79
80
81    def elevation_tilt(x, y):
82        return x*slope
83       
84    domain.set_quantity('stage', elevation_tilt)
85    domain.set_quantity('friction', friction) 
86    domain.set_quantity('elevation',elevation_tilt)
87
88    print 'Available boundary tags', domain.get_boundary_tags()
89    domain.set_region('dam','stage',inital_depth,
90                                 location = 'unique vertices') 
91
92    Br = Reflective_boundary(domain)
93    Bd = Dirichlet_boundary([0,0,0])  # to drain the water out.
94    domain.set_boundary( {'wall': Br, 'edge': Bd} )
95
96    #-------------------------------------------------------------------------
97    # Evolve system through time
98    #-------------------------------------------------------------------------
99    t0 = time.time()
100
101    for t in domain.evolve(yieldstep = 0.1, finaltime = 10): #enter timestep and final time
102        domain.write_time()
103   
104        print 'That took %.2f seconds' %(time.time()-t0)
105        print 'finished'
106       
107    points = [[0.4,0.2], [0.45,0.2]]
108    #-------------------------------------------------------------------------
109    # Calculate gauge info
110    #-------------------------------------------------------------------------
111    print "name",project.basename +".sww"
112    interpolate_sww2csv( project.basename +".sww",
113                        points,
114                        project.depth_filename ,
115                        project.velocity_x_filename,
116                        project.velocity_y_filename )
117
118
119#-------------------------------------------------------------
120if __name__ == "__main__":
121    main()
122   
Note: See TracBrowser for help on using the repository browser.