source: inundation/validation/Completed/LWRU2/lwru2.py @ 1835

Last change on this file since 1835 was 1835, checked in by ole, 19 years ago

Implemented tms file format (like sww without x, y)
Fixed broken test in sww2ers

File size: 4.6 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
13#use_variable_mesh = True #Use large variable mesh generated by create_mesh.py
14use_variable_mesh = False #Use small structured mesh
15
16
17import sys
18from os import sep
19sys.path.append('..'+sep+'..'+sep)
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 pyvolution.shallow_water import Domain, Reflective_boundary,\
79     File_boundary, Transmissive_Momentum_Set_Stage_boundary
80from pyvolution.mesh_factory import rectangular_cross
81from pyvolution.pmesh2domain import pmesh_to_domain_instance
82from Numeric import array, zeros, Float, allclose
83import project
84from caching import cache
85
86#Preparing time boundary
87prepare_timeboundary(project.boundary_filename)
88
89#Preparing points
90from pyvolution.data_manager import xya2pts
91xya2pts(project.bathymetry_filename, verbose = True,
92        z_func = lambda z: -z)
93
94
95
96#######################
97# Create Domain
98if use_variable_mesh is True:
99    print 'Creating domain from', project.mesh_filename
100
101    domain = cache(pmesh_to_domain_instance,
102                   (project.mesh_filename, Domain),
103                   dependencies = [project.mesh_filename])
104
105else:
106    print 'Creating regular mesh'
107    N = 150
108    points, vertices, boundary = rectangular_cross(N, N/5*3,
109                                                   len1=5.448, len2=3.402)
110    print 'Creating domain'
111
112    #domain = Domain(points, vertices, boundary)
113    domain = cache(Domain, (points, vertices, boundary))
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
123domain.check_integrity()
124print 'Number of triangles = ', len(domain)
125print 'The extent is ', domain.get_extent()
126
127
128
129#######################
130# Initial Conditions
131print 'Initial values'
132
133domain.set_quantity('elevation',
134                    filename=project.bathymetry_filename[:-4] + '.pts',
135                    alpha = 0.0001,
136                    verbose = True,
137                    use_cache = True)
138
139domain.set_quantity('friction', 0.0)
140domain.set_quantity('stage', 0.0)
141
142
143
144######################
145# Boundary Conditions
146#
147print 'Boundaries'
148Br = Reflective_boundary(domain)
149
150from pyvolution.util import file_function
151function = file_function(project.boundary_filename[:-4] + '.tms',
152                         domain,
153                         verbose = True)
154Bts = Transmissive_Momentum_Set_Stage_boundary(domain, function)
155
156#Set boundary conditions
157if use_variable_mesh is True:
158    domain.set_boundary({'wave': Bts, 'wall': Br})
159else:
160    domain.set_boundary({'left': Bts, 'right': Br, 'bottom': Br, 'top': Br})
161
162
163
164#######################
165# Evolve
166import time
167t0 = time.time()
168
169for t in domain.evolve(yieldstep = 0.05, finaltime = 22.5):
170    domain.write_time()
171
172print 'That took %.2f seconds' %(time.time()-t0)
Note: See TracBrowser for help on using the repository browser.