source: development/okushiri_2005/lwru2.py @ 3563

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

Moved shallow water out from the old pyvolution directory.
All tests pass, most examples and okushiri works again.

File size: 4.7 KB
Line 
1"""Validation of the AnuGA implementation of the shallow water wave equation.
2
3This script sets up LWRU2 benchmark with initial condition stated
4
5See also
6
7http://www.cee.cornell.edu/longwave/index.cfm?page=benchmark&problem=2
8
9Depth at western boundary is d = 13.5 cm
10"""
11
12
13use_variable_mesh = True #Use large variable mesh generated by create_mesh.py
14#use_variable_mesh = False #Use small structured mesh for speed
15
16
17import sys
18from os import sep
19sys.path.append('..'+sep+'..'+sep) #FIXME: Shouldn't be necessary
20
21def prepare_timeboundary(filename):
22    """Converting benchmark 2 time series to NetCDF tms file.
23    This is a 'throw-away' code taylor made for files like
24    'Benchmark_2_input.txt' from the LWRU2 benchmark
25    """
26
27    print 'Preparing time boundary from %s' %filename
28    from Numeric import array
29
30    fid = open(filename)
31
32    #Skip first line
33    line = fid.readline()
34
35    #Read remaining lines
36    lines = fid.readlines()
37    fid.close()
38
39
40    N = len(lines)
41    T = zeros(N, Float)  #Time
42    Q = zeros(N, Float)  #Values
43
44    for i, line in enumerate(lines):
45        fields = line.split()
46
47        T[i] = float(fields[0])
48        Q[i] = float(fields[1])
49
50
51    #Create tms file
52    from Scientific.IO.NetCDF import NetCDFFile
53
54    outfile = filename[:-4] + '.tms'
55    print 'Writing to', outfile
56    fid = NetCDFFile(outfile, 'w')
57
58    fid.institution = 'Geoscience Australia'
59    fid.description = 'Input wave for Benchmark 2'
60    fid.starttime = 0.0
61    fid.createDimension('number_of_timesteps', len(T))
62    fid.createVariable('time', Float, ('number_of_timesteps',))
63    fid.variables['time'][:] = T
64
65    fid.createVariable('stage', Float, ('number_of_timesteps',))
66    fid.variables['stage'][:] = Q[:]
67
68    fid.createVariable('xmomentum', Float, ('number_of_timesteps',))
69    fid.variables['xmomentum'][:] = 0.0
70
71    fid.createVariable('ymomentum', Float, ('number_of_timesteps',))
72    fid.variables['ymomentum'][:] = 0.0
73
74    fid.close()
75
76
77# Module imports
78from anuga.shallow_water import Domain, Reflective_boundary,\
79     File_boundary, Transmissive_Momentum_Set_Stage_boundary
80
81from anuga.abstract_2d_finite_volumes.mesh_factory import rectangular_cross
82from anuga.abstract_2d_finite_volumes.pmesh2domain import pmesh_to_domain_instance
83from Numeric import array, zeros, Float, allclose
84import project
85from caching import cache
86
87#Preparing time boundary
88prepare_timeboundary(project.boundary_filename)
89
90#Preparing points (Only when using original Benchmark_2_Bathymetry.txt file)
91#from anuga.abstract_2d_finite_volumes.data_manager import xya2pts
92#xya2pts(project.bathymetry_filename, verbose = True,
93#        z_func = lambda z: -z)
94
95
96
97#######################
98# Create Domain
99if use_variable_mesh is True:
100    print 'Creating domain from', project.mesh_filename
101
102    domain = Domain(project.mesh_filename, use_cache=True, verbose=True)
103
104else:
105    print 'Creating regular mesh'
106    N = 150
107    points, vertices, boundary = rectangular_cross(N, N/5*3,
108                                                   len1=5.448, len2=3.402)
109    print 'Creating domain'
110
111    #domain = Domain(points, vertices, boundary)
112    domain = cache(Domain, (points, vertices, boundary))
113
114
115
116
117import sys, os
118base = os.path.basename(sys.argv[0])
119domain.filename, _ = os.path.splitext(base)
120domain.default_order = 2
121domain.store = True    #Store for visualisation purposes
122
123#domain.check_integrity()
124print 'Number of triangles = ', len(domain)
125print 'The extent is ', domain.get_extent()
126print domain.statistics()
127
128
129
130#######################
131# Initial Conditions
132print 'Initial values'
133
134domain.set_quantity('elevation',
135                    filename = project.bathymetry_filename[:-4] + '.pts',
136                    alpha = 0.001,
137                    verbose = True,
138                    use_cache = True)
139
140domain.set_quantity('friction', 0.0)
141domain.set_quantity('stage', 0.0)
142
143
144
145######################
146# Boundary Conditions
147#
148print 'Boundaries'
149Br = Reflective_boundary(domain)
150
151from anuga.abstract_2d_finite_volumes.util import file_function
152function = file_function(project.boundary_filename[:-4] + '.tms',
153                         domain,
154                         verbose = True)
155Bts = Transmissive_Momentum_Set_Stage_boundary(domain, function)
156
157#Set boundary conditions
158if use_variable_mesh is True:
159    domain.set_boundary({'wave': Bts, 'wall': Br})
160else:
161    domain.set_boundary({'left': Bts, 'right': Br, 'bottom': Br, 'top': Br})
162
163
164
165#######################
166# Evolve
167import time
168t0 = time.time()
169
170for t in domain.evolve(yieldstep = 0.05, finaltime = 22.5):
171    domain.write_time()
172
173print 'That took %.2f seconds' %(time.time()-t0)
Note: See TracBrowser for help on using the repository browser.