source: anuga_work/development/LWRU1/lwru1.py @ 3970

Last change on this file since 3970 was 3846, checked in by ole, 18 years ago

Refactored references to domain.filename away.
Use

domain.set_name()
domain.get_name()

instead.

File size: 5.7 KB
Line 
1"""Example of shallow water wave equation.
2
3This script sets up a 2D version of the 1D LWRU1 benchmark with initial condition stated in the file benchmark_1.txt.
4
5See also
6
7http://www.cee.cornell.edu/longwave/index.cfm?page=benchmark&problem=1
8
9
10"""
11
12######################
13# Module imports
14
15from anuga.pyvolution.shallow_water import Domain, Reflective_boundary,\
16     Dirichlet_boundary,Transmissive_boundary, Constant_height, Constant_stage
17
18from anuga.pyvolution.mesh_factory import rectangular_cross
19from Numeric import array, zeros, Float, allclose
20
21
22#######################
23# Domain
24#
25
26
27print 'Creating domain'
28#Create basic mesh
29#
30#The initial condition extends 50km off shore
31#and 5,000m is allowed on shore for wetting
32#(only about 200m is expected, though)
33
34points, vertices, boundary = rectangular_cross(150, 15,
35                                         len1=55000, len2=5000,
36                                         origin = (-5000, 0.0))
37
38#points, vertices, boundary = rectangular_cross(100, 10,
39#                                         len1=55000, len2=5000,
40#                                         origin = (-5000, 0.0))
41
42
43#Create shallow water domain
44domain = Domain(points, vertices, boundary)
45
46domain.check_integrity()
47domain.default_order = 2
48
49#Output params
50domain.smooth = True
51domain.reduction = min  #Looks a lot better on top of steep slopes
52print "Number of triangles = ", len(domain)
53
54domain.visualise = False
55domain.store = True    #Store for visualisation purposes
56domain.format = 'sww'   #Native netcdf visualisation format
57
58import sys, os
59base = os.path.basename(sys.argv[0])
60
61basename, _ = os.path.splitext(base)
62domain.set_name(basename)
63
64#Set initial values
65def slope(x, y):
66    return -x/10
67
68
69class IC_x:
70    """
71    Read 1D initial condition and provide values at any x, y
72
73    File is assumed to list x values in the first column and
74    stage in the second.
75    """
76
77    def __init__(self, filename):
78
79        self.x = []
80        self.w = []
81        fid = open(filename)
82        for line in fid.readlines():
83            fields = line.split()
84            assert len(fields) == 2, '%s' %fields
85            self.x.append( float(fields[0]) )
86            self.w.append( float(fields[1]) )
87
88        #print 'X', self.x, len(self.x)
89        #print 'W', self.w, len(self.w)
90        #from pylab import plot, show
91        #plot(self.x, self.w)
92        #show()
93        #import sys; sys.exit()
94
95    def __call__(self, x, y):
96
97        w = zeros( len(x), Float )
98        for i in range(len(x)):
99            xi = x[i]
100
101
102            #Find slot
103
104            if xi < self.x[0]:
105                w[i] = self.w[0]
106            elif xi > self.x[-1]:
107                w[i] = self.w[-1]
108            else:
109                index = 0
110                while xi > self.x[index]: index += 1
111                while xi < self.x[index]: index -= 1
112
113                #print xi, index, self.x[index], self.w[index]
114
115                if xi == self.x[index]:
116                #if allclose(xi, self.x[index]):
117                    #Protect against case where x is the last value
118                    # - also works in general when x == self.x[i]
119                    ratio = 0
120                else:
121                    #x is now between index and index+1
122                    ratio = (xi - self.x[index])/\
123                            (self.x[index+1] - self.x[index])
124
125                #print xi, index, self.x[index], ratio
126
127                #Compute interpolated value
128                if ratio > 0:
129                    w[i] = self.w[index] +\
130                           ratio*(self.w[index+1] - self.w[index])
131                else:
132                    w[i] = self.w[index]
133
134        #print x, w
135        return w
136
137
138
139print 'Field values'
140domain.set_quantity('elevation', slope)
141domain.set_quantity('friction', 0.0)
142domain.set_quantity('stage', IC_x('lwru1_IC.txt'))
143
144#import sys; sys.exit()
145
146#print domain.quantities['stage'].centroid_values
147
148######################
149# Boundary conditions
150#
151print 'Boundaries'
152Br = Reflective_boundary(domain)
153Bt = Transmissive_boundary(domain)
154
155#Constant inflow
156Bd = Dirichlet_boundary([0.0, 0.0, 0.0])
157
158#Set boundary conditions
159domain.set_boundary({'left': Br, 'right': Br, 'bottom': Br, 'top': Br})
160
161
162#Evolve
163import time
164t0 = time.time()
165
166
167
168
169pt = []
170xes = []
171y = 2500
172x0 = -500
173step = 5
174for i in range(1000):
175    x = x0+i*step
176    xes.append(x)
177    pt.append( [x,y] )
178
179from pylab import *
180from anuga.pyvolution.least_squares import Interpolation
181
182
183V = domain.get_vertex_coordinates(obj=True) #Why?
184T = domain.get_triangles(obj=True)
185
186
187I = Interpolation(V,
188                  T,
189                  point_coordinates = pt,
190                  verbose = True)
191
192
193f = domain.quantities['elevation'].vertex_values.flat
194z = I.interpolate( f )
195
196print 'xxxxx'
197
198
199f = domain.quantities['stage'].vertex_values.flat
200y = I.interpolate( f )
201
202#ion()
203#plot(xes, y, '-b', xes, z, '-k', [-500, 50000], [0.0, 0.0], '-k')
204ion()
205clf()
206hold(True)
207plot(xes, y, '-b')
208plot(xes, z, '-k')
209plot([-500, 50000], [0.0, 0.0], '-k')
210set( gca(), Ylim=(-100,100) )
211set( gca(), Xlim=(-500,2000) )
212draw()
213ioff()
214
215#raw_input('go')
216for t in domain.evolve(yieldstep = 10, finaltime = 300.0):
217    domain.write_time()
218
219
220    f = domain.quantities['stage'].vertex_values.flat
221    y = I.interpolate( f )
222
223    clf()
224    hold(True)
225    plot(xes, y, '-b')
226    plot(xes, z, '-k')
227    plot([-500, 50000], [0.0, 0.0], '-k')
228    set( gca(), Ylim=(-100,100) )
229    set( gca(), Xlim=(-500,2000) )
230    draw()
231
232
233    #raw_input('go')
234
235    #print y[:], y.shape
236
237
238
239
240
241print 'That took %.2f seconds' %(time.time()-t0)
242show()
Note: See TracBrowser for help on using the repository browser.