source: inundation/examples/run_tsh.py @ 1809

Last change on this file since 1809 was 1809, checked in by duncan, 19 years ago

print georef info

File size: 4.8 KB
Line 
1"""Example of shallow water wave equation.
2
3Specific methods pertaining to the 2D shallow water equation
4are imported from shallow_water
5for use with the generic finite volume framework
6
7A example of running this program is;
8python run_tsh.py visualise hill.tsh 0.05 1
9"""
10
11######################
12# Module imports
13#
14
15import sys
16from os import sep, path
17sys.path.append('..'+sep+'pyvolution')
18
19from shallow_water import Domain, Reflective_boundary, Dirichlet_boundary,\
20     Transmissive_boundary, Time_boundary, Add_value_to_region
21from mesh_factory import rectangular
22from pmesh2domain import pmesh_to_domain_instance
23
24from Numeric import array
25
26import time
27
28#from config import default_datadir
29
30######################
31# Domain
32
33import sys
34
35
36    ######NEW
37def add_x_y(x, y):
38    return x+y
39
40    ######NEW
41
42usage = "usage: %s ['visual'|'non-visual'] pmesh_file_name yieldstep finaltime" %         path.basename(sys.argv[0])
43
44if len(sys.argv) < 4:
45    print usage
46else:
47    if sys.argv[1][0] == "n" or sys.argv[1][0] == "N":
48        visualise = False
49    else:   
50        visualise = True
51    filename = sys.argv[2]
52    yieldstep = float(sys.argv[3])
53    finaltime = float(sys.argv[4])
54   
55    print 'Creating domain from', filename
56    domain = pmesh_to_domain_instance(filename, Domain)
57    print "Number of triangles = ", len(domain)
58    print "domain.geo_reference",domain.geo_reference
59    domain.checkpoint = False #True
60    domain.default_order = 1
61    domain.visualise = visualise
62    domain.smooth = True
63    domain.set_datadir('.')
64
65    if (domain.visualise):
66        domain.store = False  #True    #Store for visualisation purposes
67    else:
68        domain.store = True  #True    #Store for visualisation purposes
69        domain.format = 'sww'   #Native netcdf visualisation format
70   
71        file_path, filename = path.split(filename)
72        filename, ext = path.splitext(filename)
73        if domain.smooth is True:
74            s = 'smooth'
75        else:
76            s = 'nonsmooth'       
77        domain.filename = filename + '_' + s + '_ys'+ str(yieldstep) + \
78                          '_ft' + str(finaltime)
79        print "Output being written to " + domain.get_datadir() + sep + \
80              domain.filename + "." + domain.format
81
82
83    #Set friction
84    manning = 0.07
85    inflow_stage = 10.0
86
87    ######NEW
88    domain.set_quantity('friction', manning)
89    domain.set_quantity('stage', add_x_y)
90    domain.set_quantity('elevation',
91                            domain.quantities['stage'].vertex_values+ \
92                            domain.quantities['elevation'].vertex_values)
93    domain.set_quantity('stage', 0.0)
94    ######NEW
95
96    ######################
97    # Boundary conditions
98    #
99    print 'Boundaries'
100    reflective = Reflective_boundary(domain)
101    Bt = Transmissive_boundary(domain)
102
103    #Constant inflow
104    Bd = Dirichlet_boundary(array([10, 0.0, 0.0]))
105
106    #Time dependent inflow
107    from math import sin, pi
108    Bw = Time_boundary(domain=domain,
109                       f=lambda x: array([(1 + sin(x*pi/4))*\
110                        (inflow_stage*(sin(2.5*x*pi)+0.7)),0,0]))
111
112
113    print 'Available boundary tags are', domain.get_boundary_tags()
114
115    #Set boundary conditions
116   
117    tags = {}
118    tags['left'] = Bw
119    tags['1'] = Time_boundary(domain=domain,
120                       f=lambda x: array([(1 + sin(x*pi/4))*\
121                        (0.15*(sin(2.5*x*pi)+0.7)),0,0]))
122    tags['wave'] = Bd
123    tags['wave'] = Time_boundary(domain=domain,
124                       f=lambda x: array([(1 + sin(x*pi/4))*\
125                        (0.15*(sin(2.5*x*pi)+0.7)),0,0]))
126    tags['internal'] = None
127    tags['levee'] = None
128    tags['0'] = reflective
129    tags['wall'] = reflective
130    tags['external'] = reflective 
131    tags['exterior'] = reflective
132    tags['open'] = Bd 
133    tags['opening'] = None 
134
135    domain.set_boundary(tags)
136
137    # region tags
138    #domain.set_region(Set_region('mound', 'elevation', 100, location='unique vertices'))
139    domain.set_region(Add_value_to_region('floodplain', 'elevation', 100, location='unique vertices', initial_quantity='elevation'))
140    domain.set_region(Add_value_to_region('mound', 'elevation', 100, location='unique vertices', initial_quantity='elevation'))
141    domain.set_region(Add_value_to_region('reservoir', 'stage', 100.0,initial_quantity='elevation'))
142   
143    #print domain.quantities['elevation'].vertex_values
144    #print domain.quantities['stage'].vertex_values
145         
146    domain.check_integrity()
147
148    ######################
149    #Evolution
150    t0 = time.time()
151    for t in domain.evolve(yieldstep = yieldstep, finaltime = finaltime):
152        domain.write_time()
153   
154    print 'That took %.2f seconds' %(time.time()-t0)
155
156   
Note: See TracBrowser for help on using the repository browser.