source: anuga_work/development/demos/netherlands.py @ 5153

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

Updated the old netherlands example in preparation for Steve's new deep-water limiter.

File size: 4.8 KB
Line 
1"""Example of shallow water wave equation.
2
3This is called Netherlands because it shows a dam with a gap in it and
4stylised housed behind it and below the water surface.
5
6"""
7
8#------------------------------------------------------------------------------
9# Import necessary modules
10#------------------------------------------------------------------------------
11
12from anuga.shallow_water import Domain
13from anuga.shallow_water import Reflective_boundary, Dirichlet_boundary
14from anuga.abstract_2d_finite_volumes.mesh_factory import rectangular_cross
15import os
16
17#from anuga.visualiser import RealtimeVisualiser
18#import rpdb
19#rpdb.set_active()
20
21
22#------------------------------------------------------------------------------
23# Setup computational domain
24#------------------------------------------------------------------------------
25
26N = 150 # size = 45000
27N = 130 # size = 33800
28N = 600 # Size = 720000
29N = 100
30
31points, elements, boundary = rectangular_cross(N, N)
32domain = Domain(points, elements, boundary, use_inscribed_circle=True)
33
34domain.check_integrity()
35
36domain.set_name(os.path.splitext(__file__)[0])
37domain.set_timestepping_method('rk3')
38domain.set_default_order(2)
39domain.set_store_vertices_uniquely(True) # Store as internally represented
40domain.tight_slope_limiters = True
41print domain.statistics()
42
43
44# Setup order and all the beta's for the limiters (these should become defaults
45
46#domain.beta_w      = 1.0
47#domain.beta_w_dry  = 0.2
48#domain.beta_uh     = 1.0
49#domain.beta_uh_dry = 0.2
50#domain.beta_vh     = 1.0
51#domain.beta_vh_dry = 0.2
52
53#domain.alpha_balance = 100.0
54
55
56
57#------------------------------------------------------------------------------
58# Setup initial conditions
59#------------------------------------------------------------------------------
60
61class Weir:
62    """Set a bathymetry for simple weir with a hole.
63    x,y are assumed to be in the unit square
64    """
65
66    def __init__(self, stage):
67        self.inflow_stage = stage
68
69    def __call__(self, x, y):
70        from Numeric import zeros, Float
71
72        N = len(x)
73        assert N == len(y)
74
75        z = zeros(N, Float)
76        for i in range(N):
77            z[i] = -x[i]/20  # General slope
78
79            # Flattish bit to the left
80            if x[i] <= 0.3:
81                #z[i] = -x[i]/5
82                z[i] = -x[i]/20
83
84
85            # Weir
86            if x[i] > 0.3 and x[i] < 0.4:
87                z[i] = -x[i]/20+1.2
88
89            # Dip
90            #if x[i] > 0.6 and x[i] < 0.9:
91            #    z[i] = -x[i]/20-0.5  #-y[i]/5
92
93            # Hole in weir
94            #if x[i] > 0.3 and x[i] < 0.4 and y[i] > 0.2 and y[i] < 0.4:
95            if x[i] > 0.3 and x[i] < 0.4 and y[i] > 0.4 and y[i] < 0.6:
96                #z[i] = -x[i]/5
97                z[i] = -x[i]/20
98
99            # Poles
100            #if x[i] > 0.65 and x[i] < 0.8 and y[i] > 0.55 and y[i] < 0.65 or\
101            #       x[i] > 0.75 and x[i] < 0.9 and y[i] > 0.35 and y[i] < 0.45:
102            #    z[i] = -x[i]/20+0.4
103
104            if (x[i] - 0.72)**2 + (y[i] - 0.6)**2 < 0.05**2:# or\
105                   #x[i] > 0.75 and x[i] < 0.9 and y[i] > 0.35 and y[i] < 0.45:
106                z[i] = -x[i]/20+0.4
107
108
109            # Wall
110            if x[i] > 0.995:
111                z[i] = -x[i]/20+0.3
112
113        return z/2
114
115
116inflow_stage = 0.5
117manning = 0.0
118
119domain.set_quantity('elevation', Weir(inflow_stage))
120domain.set_quantity('friction', manning)
121domain.set_quantity('stage', expression='elevation + 0.0')
122
123
124#------------------------------------------------------------------------------
125# Setup boundary conditions
126#------------------------------------------------------------------------------
127
128Br = Reflective_boundary(domain)
129Bd = Dirichlet_boundary([inflow_stage, 0.0, 0.0]) # Constant inflow
130domain.set_boundary({'left': Bd, 'right': Br, 'bottom': Br, 'top': Br})
131
132
133#------------------------------------------------------------------------------
134# Evolve system through time
135#------------------------------------------------------------------------------
136
137
138if N <= 150:
139    # Initialise real-time visualiser
140
141    pass
142    #vis = RealtimeVisualiser(domain)
143    #vis.render_quantity_height("elevation", dynamic=False)
144    #vis.render_quantity_height("stage", dynamic=True)
145    #vis.colour_height_quantity('stage', (0.0, 0.0, 0.8))
146    #vis.start()
147
148
149
150import time
151t0 = time.time()
152
153for t in domain.evolve(yieldstep = 0.005, finaltime = None):
154    print domain.timestepping_statistics()
155    print domain.quantities['stage'].get_values(location='centroids',
156                                                indices=[0])
157                                               
158    #vis.update()                                               
159    #time.sleep(0.1)
160    #raw_input('pause>')
161    #V.update_quantity('stage')
162    #rpdb.set_active()
163    #integral_label.text='Integral=%10.5e'%domain.quantities['stage'].get_integral()
164#vis.evolveFinished()
165
166print 'That took %.2f seconds' %(time.time()-t0)
167
168#vis.join()
Note: See TracBrowser for help on using the repository browser.