source: anuga_validation/okushiri_2005/okushiri_parallel.py @ 3645

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

Work on okushiri parallel

File size: 3.9 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# Preparing time boundary
85#-------------------------
86if myid == 0:
87    prepare_timeboundary(project.boundary_filename)
88
89
90#-------------------------
91# Create Domain
92#-------------------------
93domain = Domain(project.mesh_filename, use_cache=True, verbose=True)
94domain.set_name('okushiri')
95domain.set_default_order(2)
96print domain.statistics()
97
98
99#-------------------------
100# Initial Conditions
101#-------------------------
102domain.set_quantity('friction', 0.0)
103domain.set_quantity('stage', 0.0)
104domain.set_quantity('elevation',
105                    filename = project.bathymetry_filename[:-4] + '.pts',
106                    alpha = 0.02,                   
107                    verbose = True,
108                    use_cache = True)
109
110#-------------------------
111# Boundary Conditions
112#-------------------------
113Br = Reflective_boundary(domain)
114
115# Note this boundary cannot be fully specified here and
116# distributed without problem.
117#Bts = Transmissive_Momentum_Set_Stage_boundary(domain, function)
118
119domain.set_boundary({'wave': None,  # Bind this one later
120                     'wall': Br})
121
122
123#-------------------------
124# Distribute domain
125#-------------------------
126if numprocs > 1:
127    domain = distribute(domain)
128
129function = file_function(project.boundary_filename[:-4] + '.tms',
130                         domain, 
131                         verbose=True)
132
133Bts = Transmissive_Momentum_Set_Stage_boundary(domain, function)
134domain.modify_boundary({'wave': Bts})
135
136#-------------------------
137# Evolve through time
138#-------------------------
139import time
140t0 = time.time()
141
142for t in domain.evolve(yieldstep = 0.05, finaltime = 22.5):
143    domain.write_time()
144
145print 'That took %.2f seconds' %(time.time()-t0)
Note: See TracBrowser for help on using the repository browser.