source: anuga_work/development/demos/run_tsh.py @ 5236

Last change on this file since 5236 was 4713, checked in by steve, 18 years ago

Implemented 2nd and 3rd order timestepping via the domain.set_timestepping_method method

File size: 5.1 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 n hill.tsh 0.05 1
9"""
10
11######################
12# Module imports
13#
14
15from Numeric import array
16import time
17import sys
18from os import sep, path
19
20from anuga.shallow_water import Domain, Reflective_boundary, \
21     Dirichlet_boundary, Transmissive_boundary, Time_boundary
22from anuga.abstract_2d_finite_volumes.region import Add_value_to_region, \
23     Set_region
24from anuga.visualiser import RealtimeVisualiser
25
26
27
28#from anuga.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       
52    filename = sys.argv[2]
53    yieldstep = float(sys.argv[3])
54    finaltime = float(sys.argv[4])
55   
56    print 'Creating domain from', filename
57    domain = Domain(filename)
58
59    # check if the visualiser will work
60    try:
61        xx = RealtimeVisualiser(domain)
62    except:
63        print "Warning: Error in RealtimeVisualiser. Could not visualise."
64        visualise = False
65   
66    print "Number of triangles = ", len(domain)
67    print "domain.geo_reference",domain.geo_reference
68    domain.checkpoint = False #True
69    domain.default_order = 1
70    domain.smooth = True
71    domain.set_datadir('.')
72
73    if (visualise):
74        domain.store = False  #True    #Store for visualisation purposes
75    else:
76        domain.store = True  #True    #Store for visualisation purposes
77        domain.format = 'sww'   #Native netcdf visualisation format
78   
79        file_path, filename = path.split(filename)
80        filename, ext = path.splitext(filename)
81        if domain.smooth is True:
82            s = 'smooth'
83        else:
84            s = 'nonsmooth'       
85        domain.set_name(filename + '_' + s + '_ys'+ str(yieldstep) + \
86                        '_ft' + str(finaltime))
87        print "Output being written to " + domain.get_datadir() + sep + \
88              domain.get_name() + "." + domain.format
89
90
91    #Set friction
92    manning = 0.07
93    inflow_stage = 10.0
94
95   
96    domain.set_quantity('friction', manning)
97
98    #domain.set_quantity('stage', add_x_y)
99    #domain.set_quantity('elevation',
100    #                        domain.quantities['stage'].vertex_values+ \
101    #                        domain.quantities['elevation'].vertex_values)
102    #domain.set_quantity('stage', 0.0)
103   
104
105    ######################
106    # Boundary conditions
107    #
108    print 'Boundaries'
109    reflective = Reflective_boundary(domain)
110    Bt = Transmissive_boundary(domain)
111
112    #Constant inflow
113    Bd = Dirichlet_boundary(array([3, 0.0, 0.0]))
114
115    #Time dependent inflow
116    from math import sin, pi
117    Bw = Time_boundary(domain=domain,
118                       f=lambda x: array([(1 + sin(x*pi/4))*\
119                        (inflow_stage*(sin(2.5*x*pi)+0.7)),0,0]))
120
121
122    print 'Available boundary tags are', domain.get_boundary_tags()
123
124    #Set boundary conditions
125   
126    tags = {}
127    tags['left'] = Bw
128    tags['1'] = Bd
129   
130    tags['wave'] = Bd
131    tags['wave'] = Time_boundary(domain=domain,
132                       f=lambda x: array([(1 + sin(x*pi/4))*\
133                        (0.15*(sin(2.5*x*pi)+0.7)),0,0]))
134    tags['internal'] = None
135    tags['levee'] = None
136    tags['0'] = reflective
137    tags['wall'] = reflective
138    tags['external'] = reflective 
139    tags['exterior'] = reflective
140    tags['open'] = Bd 
141    tags['opening'] = None 
142
143    domain.set_boundary(tags)
144
145    # region tags
146   
147    domain.set_region(Set_region('slow', 'friction', 20, location='unique vertices'))
148    domain.set_region(Set_region('silo', 'elevation', 20, location='unique vertices'))
149    domain.set_region(Set_region('wet', 'elevation', 0, location='unique vertices'))
150    domain.set_region(Set_region('dry', 'elevation', 2, location='unique vertices'))
151    domain.set_region(Add_value_to_region('wet', 'stage', 1.5, location='unique vertices', initial_quantity='elevation'))
152    domain.set_region(Add_value_to_region('dry', 'stage', 0, location='unique vertices', initial_quantity='elevation'))
153   
154    #print domain.quantities['elevation'].vertex_values
155    #print domain.quantities['stage'].vertex_values
156         
157    domain.check_integrity()
158   
159    # prepare the visualiser
160    if visualise is True:
161        v = RealtimeVisualiser(domain)
162        v.render_quantity_height('elevation', dynamic=False)
163        v.render_quantity_height('stage', dynamic=True)
164        v.colour_height_quantity('stage', (0.0, 0.0, 0.8))
165        v.start()
166    ######################
167    #Evolution
168    t0 = time.time()
169    for t in domain.evolve(yieldstep = yieldstep, finaltime = finaltime):
170        domain.write_time()
171        if visualise is True:
172            v.update()
173    if visualise is True:
174        v.evolveFinished()
175
176   
177    print 'That took %.2f seconds' %(time.time()-t0)
178
179   
Note: See TracBrowser for help on using the repository browser.