source: anuga_validation/automated_validation_tests/urs_mux_files_validation/run_Bf.py @ 6160

Last change on this file since 6160 was 6160, checked in by rwilson, 15 years ago

Change Numeric imports to general form - ready to change to NumPy?.

File size: 9.2 KB
RevLine 
[6160]1"""
2
3this test compares the gauge data from a boundary sww file which is precomputed and
4stored in this directory, and the same gauge location from a simple model that uses
5the same boundary sww file for its boundary conditions
6
7Basically tests that file_boundary and evolve is working correctly
8"""
9
10#------------------------------------------------------------------------------
11# Import necessary modules
12#------------------------------------------------------------------------------
13
14# Standard modules
15from os import sep,getcwd, access, F_OK, mkdir, getenv
16from os.path import dirname, basename,abspath
17from shutil import copy
18import time, sys, os, tempfile
19
20# Related major packages
21from anuga.shallow_water import Domain,Dirichlet_boundary,File_boundary,Transmissive_boundary, Field_boundary
22import Numeric as num
23from anuga.pmesh.mesh_interface import create_mesh_from_regions
24from anuga.abstract_2d_finite_volumes.util import start_screen_catcher, copy_code_files, sww2timeseries, get_data_from_file
25# Application specific imports
26
27#------------------------------------------------------------------------------
28# Copy scripts to time stamped output directory and capture screen
29# output to file
30#------------------------------------------------------------------------------
31out_dir = dirname(abspath(__file__))+sep
32fileName="temp"
33
34meshes_dir_name = 'small.tsh' # this will be local
35
36tide = 2.4
37
38#--------------------------------------------------------------------------
39# Create the triangular mesh based on overall clipping polygon with a
40# tagged
41# boundary and interior regions defined in project.py along with
42# resolutions (maximal area of per triangle) for each polygon
43#--------------------------------------------------------------------------
44
45all = [[465184,7764500],[470397,7764510],[470407,7758988],[465195,7758979]]
46# N, W, S, E
47
48create_mesh_from_regions(all,
49                             boundary_tags={'ocean': [ 0],'side': [1, 3], 'back': [2]},
50                             maximum_triangle_area=50000,
51                             filename=meshes_dir_name,
52                             use_cache=False,
53                             )
54
55#-------------------------------------------------------------------------
56# Setup computational domain
57#-------------------------------------------------------------------------
58#print 'Setup computational domain'
59domain = Domain( meshes_dir_name, verbose=True)
60
61from anuga.shallow_water.data_manager import urs2sww
62boundaries_dir_name = 'o_test'
63
64# convert MUX urs files to an SWW file output
65urs2sww(boundaries_dir_name,boundaries_dir_name,
66           mint=9200, maxt= 11200,
67           fail_on_NaN= False,
68           verbose=True)
69   
70#-------------------------------------------------------------------------
71# Setup initial conditions
72#-------------------------------------------------------------------------
73
74#print 'Start Set quantity'
75
76domain.set_quantity('elevation', -42.3)
77
78#print 'Finished Set quantity'
79
80#------------------------------------------------------
81# Set domain parameters
82#------------------------------------------------------
83
84domain.set_quantity('stage', tide)
85domain.set_quantity('friction', 0.01)
86domain.set_name(fileName)
87domain.set_datadir(out_dir)
88
89#-------------------------------------------------------------------------
90# Setup boundary conditions
91#-------------------------------------------------------------------------
92Bf = Field_boundary(out_dir + 'o_test.sww',
93                  domain, mean_stage=tide, use_cache=True, verbose=False)
94
95#Br = Reflective_boundary(domain)
96Bt = Transmissive_boundary(domain)
97Bd = Dirichlet_boundary([tide,0,0])
98
99
100domain.set_boundary({'back': Bt,
101                     'side': Bd,
102                    'ocean': Bf
103                    })
104#----------------------------------------------------------------------------
105# Evolve system through time
106#----------------------------------------------------------------------------
107
108t0 = time.time()
109
110for t in domain.evolve(yieldstep = 60, finaltime = 1920):
111    domain.write_time()
112    domain.write_boundary_statistics()
113
114
115#----------------------------------------------------------------------------
116#Gets timeseries from boundary sww and evolved sww
117#----------------------------------------------------------------------------
118home = getenv('INUNDATIONHOME') #Sandpit's parent dir   
119data = 'data'
120state = 'western_australia'
121scenario_name = 'dampier.sww'
122
123scenario = 'dampier_tsunami_scenario_2006'
124an = 'anuga'
125bo = 'boundaries'
126
127run_time = 'blank'
128production_dirs = {run_time: 'URS evolved data'#,
129                   }
130
131topo = 'topographies'
132out = 'outputs'
133urs = 'urs'
134gridded = '1_10000'
135
136
137gauge_boundary_filename = 'gauges_time_series_b_near_top.csv'
138gauge_evolved_filename = 'gauges_time_series_near_top.csv'
139
140boundary_dir_filename = os.path.join(out_dir,gauge_boundary_filename)
141#print'boundary_dir_filename',boundary_dir_filename
142
143evolved_dir_filename= os.path.join(out_dir,gauge_evolved_filename)
144
145swwfiles = {}
146swwfile = out_dir + fileName + '.sww'
147swwfiles[swwfile] = run_time
148#----------------------------------------------------------------------------
149#get timeseries data from evolved sww file       
150#----------------------------------------------------------------------------
151texname, elev_output = sww2timeseries(swwfiles,
152                                      out_dir+sep+"gauges.txt",
153#                                      out_dir+sep+"gauges.csv",
154                                      production_dirs,
155                                      report = False,
156                                      plot_quantity = ['stage', 'xmomentum', 'ymomentum'],
157                                      surface = False,
158                                      time_min = None,
159                                      time_max = None,
160                                      title_on = False,
161                                      use_cache = True,
162                                      verbose = False)
163                                     
164swwfiles = {}
165swwfile = out_dir + boundaries_dir_name + '.sww'
166swwfiles[swwfile] = run_time
167#print"swwfiles",swwfiles,"shallow_water"
168
169#----------------------------------------------------------------------------
170#get timeseries data from boundary sww file       
171#----------------------------------------------------------------------------
172texname, elev_output = sww2timeseries(swwfiles,
173                                      out_dir+sep+"boundary_gauges.txt",
174#                                      out_dir+sep+"boundary_gauges.csv",
175                                      production_dirs,
176                                      report = False,
177                                      plot_quantity = ['stage', 'xmomentum', 'ymomentum'],
178                                      surface = False,
179                                      time_min = None,
180                                      time_max = None,
181                                      title_on = False,
182                                      use_cache = True,
183                                      verbose = False)
184
185#----------------------------------------------------------------------------
186#get_data_from file returns a e_data which is an array containing the
187#data from the file it read
188#----------------------------------------------------------------------------
189
190e_header, e_data = get_data_from_file(evolved_dir_filename)
191
192# assign columns from array to single vector(arrays)
193e_time = e_data[:,0]
194e_stage = e_data[:,1]
195e_momentum = e_data[:,2]
196e_speed = e_data[:,3]
197e_elevation = e_data[:,4]
198
199#----------------------------------------------------------------------------
200#get_data_from file returns a b_data which is an array containing the
201#data from the file it read
202#----------------------------------------------------------------------------
203
204b_header, b_data = get_data_from_file(boundary_dir_filename)
205
206# assign columns from array to single vector(arrays)
207b_time = b_data[:,0]
208b_stage = b_data[:,1]
209b_momentum = b_data[:,2]
210b_speed = b_data[:,3]
211b_elevation = b_data[:,4]
212
213#----------------------------------------------------------------------------
214#compares the 2 models, make the vector of data from the boundary sww
215#file have the same time interval as the evolved sww vector. this
216#allows them to be compared with the allclose statment at the end
217#----------------------------------------------------------------------------
218
219j=0
220k=0
221#from boundary sww
222b_sample = []
223#from evolved sww
224e_sample = []
225for i in range(len(b_time)):
226#    if j<(len(e_time)) and b_time[i] == e_time[j]:
227    if j<(len(e_time)) and k<(len(e_time)) and b_time[i] == e_time[j]:
228        b_sample.append(float(tide+b_stage[i]))
229        e_sample.append(float(e_stage[k]))
230#        if k <len(e_time): print 'time e equal b:', b_time[i],i, j, b_stage[i], b_sample[i], e_stage[k],(len(e_time)-1)
231        j = j +1
232        k = k +1
233
234
235e_stage = e_stage.tolist()
236e_stage.pop()
237
238
239print "test successful" 
240print 'fileName',fileName+'.sww'
241os.remove(meshes_dir_name)
242#----------------------------------------------------------------------------
243# Compares the two time series
244#----------------------------------------------------------------------------
245assert num.allclose (b_sample, e_sample, 0.5, 0.5)
246
247
248
249
250
251
252
253                             
254                             
255                                         
256
257
258
259
260
Note: See TracBrowser for help on using the repository browser.