source: anuga_validation/okushiri_2005/okushiri_parallel.py @ 3640

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

Added parallel version of okushiri validation.

File size: 3.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
11This a parallel version.
12"""
13
14# Module imports
15from Numeric import array, zeros, Float, allclose
16
17from anuga.shallow_water import Domain
18from anuga.shallow_water import Reflective_boundary
19from anuga.shallow_water import Transmissive_Momentum_Set_Stage_boundary
20from anuga.abstract_2d_finite_volumes.util import file_function
21
22from anuga_parallel.parallel_api import myid, numprocs, distribute
23
24import project
25
26
27
28def prepare_timeboundary(filename):
29    """Converting benchmark 2 time series to NetCDF tms file.
30    This is a 'throw-away' code taylor made for files like
31    'Benchmark_2_input.txt' from the LWRU2 benchmark
32    """
33
34    print 'Preparing time boundary from %s' %filename
35    fid = open(filename)
36
37    #Skip first line
38    line = fid.readline()
39
40    #Read remaining lines
41    lines = fid.readlines()
42    fid.close()
43
44
45    N = len(lines)
46    T = zeros(N, Float)  #Time
47    Q = zeros(N, Float)  #Values
48
49    for i, line in enumerate(lines):
50        fields = line.split()
51
52        T[i] = float(fields[0])
53        Q[i] = float(fields[1])
54
55
56    #Create tms file
57    from Scientific.IO.NetCDF import NetCDFFile
58
59    outfile = filename[:-4] + '.tms'
60    print 'Writing to', outfile
61    fid = NetCDFFile(outfile, 'w')
62
63    fid.institution = 'Geoscience Australia'
64    fid.description = 'Input wave for Benchmark 2'
65    fid.starttime = 0.0
66    fid.createDimension('number_of_timesteps', len(T))
67    fid.createVariable('time', Float, ('number_of_timesteps',))
68    fid.variables['time'][:] = T
69
70    fid.createVariable('stage', Float, ('number_of_timesteps',))
71    fid.variables['stage'][:] = Q[:]
72
73    fid.createVariable('xmomentum', Float, ('number_of_timesteps',))
74    fid.variables['xmomentum'][:] = 0.0
75
76    fid.createVariable('ymomentum', Float, ('number_of_timesteps',))
77    fid.variables['ymomentum'][:] = 0.0
78
79    fid.close()
80
81
82
83#----------------------
84
85print 'That took %.2f seconds' %(time.time()-t0)---
86# Preparing time boundary
87#-------------------------
88if myid == 0:
89    prepare_timeboundary(project.boundary_filename)
90
91
92#-------------------------
93# Create Domain
94#-------------------------
95domain = Domain(project.mesh_filename, use_cache=True, verbose=True)
96domain.set_name('okushiri')
97domain.set_default_order(2)
98print domain.statistics()
99
100
101#-------------------------
102# Initial Conditions
103#-------------------------
104domain.set_quantity('friction', 0.0)
105domain.set_quantity('stage', 0.0)
106domain.set_quantity('elevation',
107                    filename = project.bathymetry_filename[:-4] + '.pts',
108                    alpha = 0.02,                   
109                    verbose = True,
110                    use_cache = True)
111
112#-------------------------
113# Boundary Conditions
114#-------------------------
115Br = Reflective_boundary(domain)
116
117function = file_function(project.boundary_filename[:-4] + '.tms',
118                         domain, 
119                         verbose=True)
120
121# Note this boundary can be fully specified here and
122# distributed without problem.
123Bts = Transmissive_Momentum_Set_Stage_boundary(domain, function)
124
125domain.set_boundary({'wave': Bts, 'wall': Br})
126
127
128#-------------------------
129# Distribute domain
130#-------------------------
131domain = distribute(domain)
132
133
134#-------------------------
135# Evolve through time
136#-------------------------
137import time
138t0 = time.time()
139
140for t in domain.evolve(yieldstep = 0.05, finaltime = 22.5):
141    domain.write_time()
142
143print 'That took %.2f seconds' %(time.time()-t0)
Note: See TracBrowser for help on using the repository browser.