source: inundation/ga/storm_surge/validation/LWRU2/lwru2.py @ 1697

Last change on this file since 1697 was 1697, checked in by steve, 19 years ago

Found problem with File_Boundary as used in validation test LWRU2. Have setup new BC Transmissive_Momentum_Set _Stage

File size: 5.3 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
12import sys
13from os import sep
14sys.path.append('..'+sep+'..'+sep)
15
16def prepare_timeboundary(filename):
17    """Converting benchmark 2 time series to NetCDF sww file.
18    """
19
20
21    print 'Preparing time boundary from %s' %filename
22    from Numeric import array
23
24    fid = open(filename)
25
26    #Skip first line
27    line = fid.readline()
28
29    lines = fid.readlines()
30    fid.close()
31
32    N = len(lines)
33
34    T = zeros(N, Float)  #Time
35    Q = zeros(N, Float)  #Values
36
37    for i, line in enumerate(lines):
38        fields = line.split()
39
40        T[i] = float(fields[0])
41        Q[i] = float(fields[1])
42
43    #Create sww file
44    from Scientific.IO.NetCDF import NetCDFFile
45
46    outfile = filename[:-4] + '.sww'
47    print 'Writing to', outfile
48    fid = NetCDFFile(outfile, 'w')
49
50    fid.institution = 'Geoscience Australia'
51    fid.description = 'Input wave for Benchmark 2'
52    fid.starttime = 0.0
53    fid.createDimension('number_of_timesteps', len(T))
54    fid.createVariable('time', Float, ('number_of_timesteps',))
55    fid.variables['time'][:] = T
56
57    fid.createVariable('stage', Float, ('number_of_timesteps',))
58    fid.variables['stage'][:] = Q[:]
59
60    fid.createVariable('xmomentum', Float, ('number_of_timesteps',))
61    #fid.variables['xmomentum'][:] = 0.005
62    #fid.variables['xmomentum'][:-1] = Q[1:]-Q[:-1]/0.05 #Too high
63    #fid.variables['xmomentum'][:-1] = Q[1:]-Q[:-1]/0.5   #Too low
64    #fid.variables['xmomentum'][:-1] = Q[1:]-Q[:-1]/0.1   #Too high
65    #fid.variables['xmomentum'][:-1] = Q[1:]-Q[:-1]/0.25   #Too high
66    #fid.variables['xmomentum'][:-1] = Q[1:]-Q[:-1]/0.35 #Too low
67    #fid.variables['xmomentum'][:-1] = Q[1:]-Q[:-1]/0.3  #Too low
68    #fid.variables['xmomentum'][-1] = 0.0
69
70
71    #fid.variables['xmomentum'][:] = 0.003
72    fid.variables['xmomentum'][:] = 0.0
73
74    fid.createVariable('ymomentum', Float, ('number_of_timesteps',))
75    fid.variables['ymomentum'][:] = 0.0
76
77    fid.close()
78
79
80# Module imports
81from pyvolution.shallow_water import Domain, Reflective_boundary,\
82            File_boundary, Transmissive_Momentum_Set_Stage_boundary
83from pyvolution.mesh_factory import rectangular
84from Numeric import array, zeros, Float, allclose
85
86picklefile = 'domain.pck'
87from cPickle import dump, load
88
89
90boundary_file_name = 'Benchmark_2_input.txt'
91bathymetry_file_name = 'Benchmark_2_Bathymetry.txt'
92N = 100
93
94try:
95    #raise Exception
96    fid = open(picklefile)
97    print 'Read pickled domain'
98    domain = load(fid)
99
100except:
101
102    #Preparing time boundary
103    prepare_timeboundary(boundary_file_name)
104
105    #Preparing points
106    from pyvolution.data_manager import xya2pts
107    xya2pts(bathymetry_file_name, verbose = True,
108            z_func = lambda z: -z)
109
110
111    #######################
112    # Domain
113    print 'Creating domain'
114    points, vertices, boundary = rectangular(N, N/5*3,
115                                             len1=5.448, len2=3.402)
116    domain = Domain(points, vertices, boundary)
117
118    domain.check_integrity()
119    domain.default_order = 2
120
121    print "Number of triangles = ", len(domain)
122    domain.store = True    #Store for visualisation purposes
123
124    import sys, os
125    base = os.path.basename(sys.argv[0])
126    domain.filename, _ = os.path.splitext(base)
127
128    #LS code to be included in set_quantity
129    from pyvolution import util, least_squares
130    points, attributes = util.read_xya(bathymetry_file_name[:-4] + '.pts')
131
132    #Fit attributes to mesh
133    vertex_attributes = least_squares.fit_to_mesh(domain.coordinates,
134                                                  domain.triangles,
135                                                  points,
136                                                  attributes['elevation'],
137                                                  alpha = 0,
138                                                  verbose=True)
139
140    print 'Initial values'
141    domain.set_quantity('elevation', vertex_attributes)
142    #domain.set_quantity('elevation', bathymetry_file_name[:-4] + '.pts')
143    domain.set_quantity('friction', 0.0)
144    domain.set_quantity('stage', 0.0)
145
146
147    fid = open('domain.pck', 'w')
148    dump(domain, fid)
149    fid.close()
150
151
152######################
153# Boundary conditions
154#
155print 'Boundaries'
156Br = Reflective_boundary(domain)
157#Bf = File_boundary(boundary_file_name[:-4] + '.sww', domain, verbose=True)
158
159#####################
160# Error with File Boundary. We need to only set
161# stage and let the scheme set the momentum
162from pyvolution.util import file_function
163function = file_function(boundary_file_name[:-4] + '.sww', domain,verbose = True)
164Bts = Transmissive_Momentum_Set_Stage_boundary(domain,function)
165
166#Set boundary conditions
167domain.set_boundary({'left': Bts, 'right': Br, 'bottom': Br, 'top': Br})
168
169#Evolve
170import time
171t0 = time.time()
172domain.visualise = True
173
174for t in domain.evolve(yieldstep = 0.05, finaltime = 22.5):
175    domain.write_time()
176    Stage_cv = domain.quantities['stage'].centroid_values
177    print Stage_cv[0]
178    print 'Function value',function(t)
179
180print 'That took %.2f seconds' %(time.time()-t0)
Note: See TracBrowser for help on using the repository browser.