source: branches/numpy_anuga_validation/automated_validation_tests/urs_mux_files_validation/run_Bf.py @ 6473

"""

this test compares the gauge data from a boundary sww file which is precomputed and 
stored in this directory, and the same gauge location from a simple model that uses 
the same boundary sww file for its boundary conditions

Basically tests that file_boundary and evolve is working correctly
"""

#------------------------------------------------------------------------------
# Import necessary modules
#------------------------------------------------------------------------------

# Standard modules
from os import sep,getcwd, access, F_OK, mkdir, getenv
from os.path import dirname, basename,abspath 
from shutil import copy
import time, sys, os, tempfile

# Related major packages
from anuga.shallow_water import Domain,Dirichlet_boundary,File_boundary,Transmissive_boundary, Field_boundary
import Numeric as num
from anuga.pmesh.mesh_interface import create_mesh_from_regions
from anuga.abstract_2d_finite_volumes.util import start_screen_catcher, copy_code_files, sww2timeseries, get_data_from_file
# Application specific imports

#------------------------------------------------------------------------------
# Copy scripts to time stamped output directory and capture screen
# output to file
#------------------------------------------------------------------------------
out_dir = dirname(abspath(__file__))+sep
fileName="temp"

meshes_dir_name = 'small.tsh' # this will be local

tide = 2.4

#--------------------------------------------------------------------------
# Create the triangular mesh based on overall clipping polygon with a
# tagged
# boundary and interior regions defined in project.py along with
# resolutions (maximal area of per triangle) for each polygon
#--------------------------------------------------------------------------

all = [[465184,7764500],[470397,7764510],[470407,7758988],[465195,7758979]]
# N, W, S, E

create_mesh_from_regions(all,
                             boundary_tags={'ocean': [ 0],'side': [1, 3], 'back': [2]},
                             maximum_triangle_area=50000,
                             filename=meshes_dir_name,
                             use_cache=False,
                             )

#-------------------------------------------------------------------------
# Setup computational domain
#-------------------------------------------------------------------------
#print 'Setup computational domain'
domain = Domain( meshes_dir_name, verbose=True)

from anuga.shallow_water.data_manager import urs2sww
boundaries_dir_name = 'o_test'

# convert MUX urs files to an SWW file output
urs2sww(boundaries_dir_name,boundaries_dir_name,
           mint=9200, maxt= 11200,
           fail_on_NaN= False,
           verbose=True)
    
#-------------------------------------------------------------------------
# Setup initial conditions
#-------------------------------------------------------------------------

#print 'Start Set quantity'

domain.set_quantity('elevation', -42.3)

#print 'Finished Set quantity'

#------------------------------------------------------
# Set domain parameters
#------------------------------------------------------ 

domain.set_quantity('stage', tide)
domain.set_quantity('friction', 0.01)
domain.set_name(fileName)
domain.set_datadir(out_dir)

#-------------------------------------------------------------------------
# Setup boundary conditions 
#-------------------------------------------------------------------------
Bf = Field_boundary(out_dir + 'o_test.sww',
                  domain, mean_stage=tide, use_cache=True, verbose=False)

#Br = Reflective_boundary(domain)
Bt = Transmissive_boundary(domain)
Bd = Dirichlet_boundary([tide,0,0])


domain.set_boundary({'back': Bt,
                     'side': Bd,
                    'ocean': Bf
                    })
#----------------------------------------------------------------------------
# Evolve system through time
#----------------------------------------------------------------------------

t0 = time.time()

for t in domain.evolve(yieldstep = 60, finaltime = 1920):
    domain.write_time()
    domain.write_boundary_statistics()


#----------------------------------------------------------------------------
#Gets timeseries from boundary sww and evolved sww
#----------------------------------------------------------------------------
home = getenv('INUNDATIONHOME') #Sandpit's parent dir   
data = 'data'
state = 'western_australia'
scenario_name = 'dampier.sww'

scenario = 'dampier_tsunami_scenario_2006'
an = 'anuga'
bo = 'boundaries'

run_time = 'blank'
production_dirs = {run_time: 'URS evolved data'#,
                   }

topo = 'topographies'
out = 'outputs'
urs = 'urs'
gridded = '1_10000'


gauge_boundary_filename = 'gauges_time_series_b_near_top.csv'
gauge_evolved_filename = 'gauges_time_series_near_top.csv'

boundary_dir_filename = os.path.join(out_dir,gauge_boundary_filename)
#print'boundary_dir_filename',boundary_dir_filename

evolved_dir_filename= os.path.join(out_dir,gauge_evolved_filename)

swwfiles = {}
swwfile = out_dir + fileName + '.sww'
swwfiles[swwfile] = run_time
#----------------------------------------------------------------------------
#get timeseries data from evolved sww file        
#----------------------------------------------------------------------------
texname, elev_output = sww2timeseries(swwfiles,
                                      out_dir+sep+"gauges.txt",
#                                      out_dir+sep+"gauges.csv",
                                      production_dirs,
                                      report = False,
                                      plot_quantity = ['stage', 'xmomentum', 'ymomentum'],
                                      surface = False,
                                      time_min = None,
                                      time_max = None,
                                      title_on = False,
                                      use_cache = True,
                                      verbose = False)
                                      
swwfiles = {}
swwfile = out_dir + boundaries_dir_name + '.sww'
swwfiles[swwfile] = run_time
#print"swwfiles",swwfiles,"shallow_water"

#----------------------------------------------------------------------------
#get timeseries data from boundary sww file        
#----------------------------------------------------------------------------
texname, elev_output = sww2timeseries(swwfiles,
                                      out_dir+sep+"boundary_gauges.txt",
#                                      out_dir+sep+"boundary_gauges.csv",
                                      production_dirs,
                                      report = False,
                                      plot_quantity = ['stage', 'xmomentum', 'ymomentum'],
                                      surface = False,
                                      time_min = None,
                                      time_max = None,
                                      title_on = False,
                                      use_cache = True,
                                      verbose = False)

#----------------------------------------------------------------------------
#get_data_from file returns a e_data which is an array containing the
#data from the file it read
#----------------------------------------------------------------------------

e_header, e_data = get_data_from_file(evolved_dir_filename)

# assign columns from array to single vector(arrays)
e_time = e_data[:,0]
e_stage = e_data[:,1]
e_momentum = e_data[:,2]
e_speed = e_data[:,3]
e_elevation = e_data[:,4]

#----------------------------------------------------------------------------
#get_data_from file returns a b_data which is an array containing the
#data from the file it read
#----------------------------------------------------------------------------

b_header, b_data = get_data_from_file(boundary_dir_filename)

# assign columns from array to single vector(arrays)
b_time = b_data[:,0]
b_stage = b_data[:,1]
b_momentum = b_data[:,2]
b_speed = b_data[:,3]
b_elevation = b_data[:,4]

#----------------------------------------------------------------------------
#compares the 2 models, make the vector of data from the boundary sww
#file have the same time interval as the evolved sww vector. this
#allows them to be compared with the allclose statment at the end
#----------------------------------------------------------------------------

j=0
k=0
#from boundary sww 
b_sample = []
#from evolved sww
e_sample = []
for i in range(len(b_time)):
#    if j<(len(e_time)) and b_time[i] == e_time[j]:
    if j<(len(e_time)) and k<(len(e_time)) and b_time[i] == e_time[j]:
        b_sample.append(float(tide+b_stage[i]))
        e_sample.append(float(e_stage[k]))
#        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)
        j = j +1
        k = k +1


e_stage = e_stage.tolist()
e_stage.pop()


print "test successful" 
print 'fileName',fileName+'.sww'
os.remove(meshes_dir_name)
#----------------------------------------------------------------------------
# Compares the two time series
#----------------------------------------------------------------------------
assert num.allclose (b_sample, e_sample, 0.5, 0.5)