source: inundation/validation/LWRU2/lwru2.py @ 1753

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

Embedded caching functionality within quantity.set_values and modified validation example lwru2.py to illustrate the advantages that can be gained from supervised caching.

File size: 4.5 KB
Line 
1"""Example of 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#read_mesh = True  #Use large unstructured mesh generated from create_mesh.py
14read_mesh = False  #Use small structured mesh
15
16import sys
17from os import sep
18sys.path.append('..'+sep+'..'+sep)
19
20def prepare_timeboundary(filename):
21    """Converting benchmark 2 time series to NetCDF sww file.
22    """
23
24
25    print 'Preparing time boundary from %s' %filename
26    from Numeric import array
27
28    fid = open(filename)
29
30    #Skip first line
31    line = fid.readline()
32
33    lines = fid.readlines()
34    fid.close()
35
36    N = len(lines)
37
38    T = zeros(N, Float)  #Time
39    Q = zeros(N, Float)  #Values
40
41    for i, line in enumerate(lines):
42        fields = line.split()
43
44        T[i] = float(fields[0])
45        Q[i] = float(fields[1])
46
47    #Create sww file
48    from Scientific.IO.NetCDF import NetCDFFile
49
50    outfile = filename[:-4] + '.sww'
51    print 'Writing to', outfile
52    fid = NetCDFFile(outfile, 'w')
53
54    fid.institution = 'Geoscience Australia'
55    fid.description = 'Input wave for Benchmark 2'
56    fid.starttime = 0.0
57    fid.createDimension('number_of_timesteps', len(T))
58    fid.createVariable('time', Float, ('number_of_timesteps',))
59    fid.variables['time'][:] = T
60
61    fid.createVariable('stage', Float, ('number_of_timesteps',))
62    fid.variables['stage'][:] = Q[:]
63
64    fid.createVariable('xmomentum', Float, ('number_of_timesteps',))
65    fid.variables['xmomentum'][:] = 0.0
66
67    fid.createVariable('ymomentum', Float, ('number_of_timesteps',))
68    fid.variables['ymomentum'][:] = 0.0
69
70    fid.close()
71
72
73# Module imports
74from pyvolution.shallow_water import Domain, Reflective_boundary,\
75            File_boundary, Transmissive_Momentum_Set_Stage_boundary
76from pyvolution.mesh_factory import rectangular
77from pyvolution.pmesh2domain import pmesh_to_domain_instance
78from Numeric import array, zeros, Float, allclose
79import filenames
80from caching import cache
81
82#Preparing time boundary
83prepare_timeboundary(filenames.boundary_filename)
84
85#Preparing points
86from pyvolution.data_manager import xya2pts
87xya2pts(filenames.bathymetry_filename, verbose = True,
88        z_func = lambda z: -z)
89
90
91#######################
92# Domain
93   
94if read_mesh is True:
95    print 'Creating domain from', filenames.mesh_filename
96    #domain = pmesh_to_domain_instance(filenames.mesh_filename, Domain)
97
98    domain = cache(pmesh_to_domain_instance,
99                   (filenames.mesh_filename, Domain),
100                   dependencies = [filenames.mesh_filename])
101
102
103   
104else:   
105    print 'Creating regular mesh'
106    N = 100
107    points, vertices, boundary = rectangular(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
116print "Number of triangles = ", len(domain)
117print 'The extent is ', domain.get_extent()
118
119
120#domain.check_integrity()
121domain.default_order = 2
122
123print "Number of triangles = ", len(domain)
124domain.store = True    #Store for visualisation purposes
125
126import sys, os
127base = os.path.basename(sys.argv[0])
128domain.filename, _ = os.path.splitext(base)
129
130print 'Initial values'
131
132domain.set_quantity('elevation',
133                    filename=filenames.bathymetry_filename[:-4] + '.pts',
134                    alpha = 0.0001,
135                    verbose = True,
136                    use_cache = True)
137                   
138domain.set_quantity('friction', 0.0)
139domain.set_quantity('stage', 0.0)
140
141
142
143######################
144# Boundary conditions
145#
146print 'Boundaries'
147Br = Reflective_boundary(domain)
148
149#####################
150# Error with File Boundary. We need to only set
151# stage and let the scheme set the momentum
152from pyvolution.util import file_function
153function = file_function(filenames.boundary_filename[:-4] + '.sww', domain,
154                         verbose = True)
155Bts = Transmissive_Momentum_Set_Stage_boundary(domain, function)
156
157#Set boundary conditions
158if read_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#Evolve
165import time
166t0 = time.time()
167#domain.visualise = True
168
169for t in domain.evolve(yieldstep = 0.05, finaltime = 22.5):
170    domain.write_time()
171    #Stage_cv = domain.quantities['stage'].centroid_values
172    #print Stage_cv[0]
173    #print 'Function value',function(t)
174
175print 'That took %.2f seconds' %(time.time()-t0)
Note: See TracBrowser for help on using the repository browser.