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

Last change on this file since 7750 was 7750, checked in by hudson, 14 years ago

Updated validation suite to work with new API - confirmed that Patong passes.

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