source: anuga_validation/okushiri_2005/lwru2.py @ 3637

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

Okushiri mesh from 18 August 2005 recreated.

File size: 4.8 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
17#import sys
18#from os import sep
19#sys.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.01,
137                    alpha = 0.02,                   
138                    verbose = True,
139                    use_cache = True)
140
141domain.set_quantity('friction', 0.0)
142domain.set_quantity('stage', 0.0)
143
144
145
146######################
147# Boundary Conditions
148#
149print 'Boundaries'
150Br = Reflective_boundary(domain)
151
152from anuga.abstract_2d_finite_volumes.util import file_function
153function = file_function(project.boundary_filename[:-4] + '.tms',
154                         domain,
155                         verbose = True)
156Bts = Transmissive_Momentum_Set_Stage_boundary(domain, function)
157
158#Set boundary conditions
159if use_variable_mesh is True:
160    domain.set_boundary({'wave': Bts, 'wall': Br})
161else:
162    domain.set_boundary({'left': Bts, 'right': Br, 'bottom': Br, 'top': Br})
163
164
165
166#######################
167# Evolve
168import time
169t0 = time.time()
170
171for t in domain.evolve(yieldstep = 0.05, finaltime = 22.5):
172    domain.write_time()
173
174print 'That took %.2f seconds' %(time.time()-t0)
Note: See TracBrowser for help on using the repository browser.