source: anuga_validation/automated_validation_tests/okushiri_tank_validation/validate_okushiri.py @ 3722

"""Automated validation of ANUGA.

This script runs the benchmark on a regular grid and verifies that
the predicted timeseries at selected gauges are close enought to the observed
timeseries provided from the experiment.

See anuga_validation/okushiri_2005 for scripts that produce simulations at
greater resolution.

"""

# Module imports
from Numeric import array, zeros, Float, allclose

from anuga.utilities.numerical_tools import ensure_numeric
from anuga.shallow_water import Domain
from anuga.shallow_water import Reflective_boundary
from anuga.shallow_water import Transmissive_Momentum_Set_Stage_boundary
from anuga.abstract_2d_finite_volumes.mesh_factory import rectangular_cross
from anuga.abstract_2d_finite_volumes.util import file_function

import project

finaltime = 22.5
timestep = 0.05

#-------------------------
# Create Domain
#-------------------------
#N = 50
N = 10
points, vertices, boundary = rectangular_cross(N, N/5*3,
                                               len1=5.448, len2=3.402)
domain = Domain(points, vertices, boundary)

domain.set_name('okushiri_automated_validation')
domain.set_default_order(2)
domain.set_minimum_storable_height(0.001)
domain.check_integrity()


#-------------------------
# Initial Conditions
#-------------------------
domain.set_quantity('friction', 0.0)
domain.set_quantity('stage', 0.0)
domain.set_quantity('elevation',
                    filename=project.bathymetry_filename[:-4] + '.pts',
                    alpha=0.02,                    
                    verbose=True,
                    use_cache=True)

#-------------------------
# Boundary Conditions
#-------------------------
Br = Reflective_boundary(domain)

function = file_function(project.boundary_filename[:-4] + '.tms',
                         domain,
                         verbose=True)
Bts = Transmissive_Momentum_Set_Stage_boundary(domain, function)
domain.set_boundary({'left': Bts, 'right': Br, 'bottom': Br, 'top': Br})


#-------------------------
# Set up validation data
#-------------------------
gauge_locations = [[0.000, 1.696]]  #Boundary
gauge_locations += [[4.521, 1.196], [4.521, 1.696], [4.521, 2.196]] #Ch 5-7-9

#Boundary input wave
filename = project.boundary_filename[:-4] + '.tms'
print 'Reading', filename
from Scientific.IO.NetCDF import NetCDFFile
from Numeric import allclose
fid = NetCDFFile(filename, 'r')#

input_time = fid.variables['time'][:]
input_stage = fid.variables['stage'][:]


# reference gauge timeseries
reference_time = []
ch5 = []
ch7 = []
ch9 = []
fid = open('output_ch5-7-9.txt')
lines = fid.readlines()
fid.close()
#for i, line in enumerate(lines[1:]):
for line in lines[1:]:
    fields = line.split()

    reference_time.append(float(fields[0]))
    ch5.append(float(fields[1])/100)   #cm2m
    ch7.append(float(fields[2])/100)   #cm2m
    ch9.append(float(fields[3])/100)   #cm2m

    if fields[0] == str(finaltime): break    
    

assert reference_time[0] == 0.0
assert reference_time[-1] == finaltime

reference_gauge_values = [ensure_numeric(input_stage),
                          ensure_numeric(ch5),
                          ensure_numeric(ch7),
                          ensure_numeric(ch9)] 



predicted_gauge_values = [[], [],[],[]]


#-------------------------
# Evolve through time 
#-------------------------
import time
t0 = time.time()

for j, t in enumerate(domain.evolve(yieldstep = timestep,
                                    finaltime = finaltime)):
    domain.write_time()

    # Record gauge values
    stage = domain.get_quantity('stage')
    for i, (x, y) in enumerate(gauge_locations):
        w = stage.get_values(interpolation_points=[[x,y]])
        predicted_gauge_values[i].append(w)

        print 'difference', x, y, w - reference_gauge_values[i][j]
    

#-------------------------
# Validate
#-------------------------
for i, (x, y) in enumerate(gauge_locations):        
    predicted_gauge_values[i] = ensure_numeric(predicted_gauge_values[i])

    
    sqsum = sum((predicted_gauge_values[i] - reference_gauge_values)**2)
    denom = sum(reference_gauge_values**2)

    print i, sqsum, sqsum/denom

    msg 'Rms error is %f' %(sqsum/denom)
    assert sqsum/denom < 0.01, msg

    try:
        from pylab import *
    except:
        pass
    else:
        ion()
        hold(False)
        plot(reference_time, reference_gauge_values[i], 'r.-',
             reference_time, predicted_gauge_values[i], 'k-')
        title('Gauge %s' %name)
        xlabel('time(s)')
        ylabel('stage (m)')    
        legend(('Observed', 'Modelled'), shadow=True, loc='upper left')
        savefig(name, dpi = 300)
        
        raw_input('Next')
        

try:
    from pylab import *
except:
    pass
else:
    show()
    

print 'That took %.2f seconds' %(time.time()-t0)