# source:anuga_work/development/Rudy_vandrie_2007/run_hydro_example.py@4440

Last change on this file since 4440 was 4440, checked in by ole, 16 years ago

Examples prepared for Rudy

File size: 3.9 KB
Line
1"""Simple water flow example using ANUGA
2
3Water flowing down a channel with more complex topography
4"""
5
6#------------------------------------------------------------------------------
7# Import necessary modules
8#------------------------------------------------------------------------------
9from anuga.abstract_2d_finite_volumes.mesh_factory import rectangular_cross
10from anuga.abstract_2d_finite_volumes.util import file_function
11from anuga.shallow_water import Domain
12from anuga.shallow_water import Reflective_boundary
13from anuga.shallow_water import Dirichlet_boundary
14from anuga.shallow_water import Time_boundary
15from anuga.shallow_water import Field_boundary
16from anuga.shallow_water.shallow_water_domain import Inflow
17
18
19#------------------------------------------------------------------------------
20# Setup computational domain
21#------------------------------------------------------------------------------
22length = 40.
23width = 5.
24dx = dy = .1           # Resolution: Length of subdivisions on both axes
25#dx = dy = .1           # Resolution: Length of subdivisions on both axes
26
27points, vertices, boundary = rectangular_cross(int(length/dx), int(width/dy),
28                                               len1=length, len2=width)
29domain = Domain(points, vertices, boundary)
30domain.set_name('hydro_example')                  # Output name
31domain.set_minimum_storable_height(0.0001)
32domain.beta_h = 0.0
33domain.limit2007 = 1
34
35
36#------------------------------------------------------------------------------
37# Setup initial conditions
38#------------------------------------------------------------------------------
39def topography(x,y):
40    """Complex topography defined by a function of vectors x and y
41    """
42
43    z = -x/10
44    #z = x*0.0
45
46    N = len(x)
47    for i in range(N):
48
49        #Step
50        if 10 < x[i] < 12:
51            z[i] += 0.4 - 0.05*y[i]
52
53        #Constriction
54        if 27 < x[i] < 29 and y[i] > 3:
55            z[i] += 2
56
57        # Pole
58        if (x[i] - 34)**2 + (y[i] - 2)**2 < 0.4**2:
59            z[i] += 2
60
61    return z
62
63
64domain.set_quantity('elevation', topography)  # Use function for elevation
65domain.set_quantity('friction', 0.01)         # Constant friction
66domain.set_quantity('stage',
67                    expression='elevation')   # Dry initial condition
68
69
70#------------------------------------------------------------------------------
71# Setup specialised forcing terms
72#------------------------------------------------------------------------------
73
75#            flow=lambda t: min(0.01*t, 0.0142)) # Tap turning up
76
78             flow=file_function('Island_Pt_Rd_Meta.tms'))
79
80domain.forcing_terms.append(hydrograph)
81
82
83#------------------------------------------------------------------------------
84# Setup boundary conditions
85#------------------------------------------------------------------------------
86
87#Hydrograph = Field_boundary('Island_Pt_Rd_Meta.tms', domain,
88#                            mean_stage=0.01)
89
90Bi = Dirichlet_boundary([0.4, 0, 0])          # Inflow
91Br = Reflective_boundary(domain)              # Solid reflective wall
92Bo = Dirichlet_boundary([-5, 0, 0])           # Outflow
93
94domain.set_boundary({'left': Br,
95                     'right': Br,
96                     'top': Br,
97                     'bottom': Br})
98
99
100#------------------------------------------------------------------------------
101# Evolve system through time
102#------------------------------------------------------------------------------
103domain.start_time = 2800
104for t in domain.evolve(yieldstep = 1, finaltime = 32700):
105    domain.write_time()
106
107
108    #if domain.get_quantity('stage').\
109    #       get_values(interpolation_points=[[10, 2.5]]) > 0:
110    #    print 'Stage > 0: Changing to outflow boundary'
111    #    domain.set_boundary({'right': Bo})
Note: See TracBrowser for help on using the repository browser.