source: anuga_validation/automated_validation_tests/UQ_runup_2006/run_dam.py @ 4751

Last change on this file since 4751 was 4751, checked in by duncan, 16 years ago

automating validation

File size: 4.7 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.region import Set_region
22from anuga.fit_interpolate.interpolate import interpolate_sww2csv, \
23     file_function
24from anuga.abstract_2d_finite_volumes.util import start_screen_catcher, \
25     copy_code_files
26from anuga.shallow_water.data_manager import csv2dict
27from numerical_tools import  err   # norm, corr, 
28
29# Application specific imports
30import create_mesh
31import project
32def main():
33     
34    slopes = [0]
35    frictions = [0.01] #can enter multiple frictions eg [0.00, 0.01]
36    inital_depths = [0.2] #can enter multiple initial depths eg [0.1, 0.2]
37    gate_positions = [0.75]
38   
39    scenario(0, 0.01, 0.02, 0.75)
40           
41
42def scenario(slope, friction, inital_depth, gate_position):
43
44    #-------------------------------------------------------------------------
45    # Create the triangular mesh
46    #-------------------------------------------------------------------------
47
48    create_mesh.generate(project.mesh_filename,
49                         gate_position,
50                         is_course=True) # this creates the mesh
51                         #is_course=False) # this creates the mesh
52
53    head,tail = path.split(project.mesh_filename)
54    #-------------------------------------------------------------------------
55    # Setup computational domain
56    #-------------------------------------------------------------------------
57    domain = Domain(project.mesh_filename, use_cache = False, verbose = True)
58   
59
60    print 'Number of triangles = ', len(domain)
61    print 'The extent is ', domain.get_extent()
62    print domain.statistics()
63
64   
65    domain.set_name(project.basename)
66    domain.set_datadir('.')
67    domain.set_quantities_to_be_stored(['stage', 'xmomentum', 'ymomentum'])
68    domain.set_minimum_storable_height(0.001)
69    domain.set_store_vertices_uniquely(True)  # for writting to sww
70
71    #-------------------------------------------------------------------------
72    # Setup initial conditions
73    #-------------------------------------------------------------------------
74
75
76    def elevation_tilt(x, y):
77        return x*slope
78       
79    domain.set_quantity('stage', elevation_tilt)
80    domain.set_quantity('friction', friction) 
81    domain.set_quantity('elevation',elevation_tilt)
82
83    print 'Available boundary tags', domain.get_boundary_tags()
84    domain.set_region('dam','stage',inital_depth,
85                                 location = 'unique vertices') 
86
87    Br = Reflective_boundary(domain)
88    Bd = Dirichlet_boundary([0,0,0])  # to drain the water out.
89    domain.set_boundary( {'wall': Br, 'edge': Bd} )
90
91    #-------------------------------------------------------------------------
92    # Evolve system through time
93    #-------------------------------------------------------------------------
94    t0 = time.time()
95
96    for t in domain.evolve(yieldstep = 0.1, finaltime = 10): #enter timestep and final time
97        domain.write_time()
98   
99        print 'That took %.2f seconds' %(time.time()-t0)
100        print 'finished'
101
102    # Load actual experimental results
103    actual,title_index_dic = csv2dict(project.actual_filename)
104   
105    gauge_locations = [[0.4,0.2]]
106    quantities = ['stage', 'elevation']
107    file_instance = file_function( project.basename +".sww",
108                                      quantities = quantities,
109                                      interpolation_points = gauge_locations,
110                                      verbose = True,
111                                      use_cache = False)
112    # create a list of the simulated_depths at the actual data times.
113    simulated_depths = []
114    for atime in actual['time']:
115        quantities_slice = file_instance(float(atime),
116                                         point_id=0)
117        depth = quantities_slice[0] - quantities_slice[1]
118        simulated_depths.append(depth)
119    flume_depths = actual["0.4:0.2"]
120    flume_depths = [float(i) for i in flume_depths]
121    # calc the norm
122    #print "map(None, simulated_depths, flume_depths)", \
123    #      map(None, simulated_depths, flume_depths)
124    normmax = 0
125    norm = err(simulated_depths,
126               flume_depths, 2, relative = True)  # 2nd norm (rel. RMS
127    if norm > normmax:
128        normmax = norm
129    print norm
130    print "norm", norm
131#-------------------------------------------------------------
132if __name__ == "__main__":
133    main()
134   
Note: See TracBrowser for help on using the repository browser.