source: anuga_validation/okushiri_2005/extract_timeseries.py @ 3895

"""Verify that simulation produced by ANUGA compares to published
validation timeseries ch5, ch7 and ch9 as well as the boundary timeseries.

RMS norm is printed and plots are produced.
"""

from Numeric import allclose
from Scientific.IO.NetCDF import NetCDFFile

from anuga.abstract_2d_finite_volumes.util import file_function
from anuga.utilities.numerical_tools import\
     ensure_numeric, cov, get_machine_precision

import project

try:
    from pylab import ion, hold, plot, title, legend, xlabel, ylabel, savefig
except:
    plotting = False
else:
    plotting = True

#plotting = False

#-------------------------
# Basic data
#-------------------------

finaltime = 22.5
timestep = 0.05

gauge_locations = [[0.000, 1.696]] # Boundary gauge
gauge_locations += [[4.521, 1.196],  [4.521, 1.696],  [4.521, 2.196]] #Ch 5-7-9
gauge_names = ['Boundary', 'ch5', 'ch7', 'ch9']

validation_data = {}
for key in gauge_names:
    validation_data[key] = []

    
#expected_covariances = {'Boundary': 5.288392008865989e-05,
#                        'ch5': 1.166748190444681e-04,
#                        'ch7': 1.121816242516758e-04,
#                        'ch9': 1.249543278366778e-04}

# old limiters
#expected_covariances = {'Boundary': 5.288601162783020386e-05,
#                        'ch5': 1.167001054284431472e-04,
#                        'ch7': 1.121474766904651861e-04,
#                        'ch9': 1.249244820847215335e-04}
#
#expected_differences = {'Boundary': 8.361144081847830638e-04,
#                        'ch5': 3.423673831653336816e-03,
#                        'ch7': 2.799962153549145211e-03,
#                        'ch9': 3.198560464876740433e-03}


#expected_covariances = {'Boundary': 5.288392008865989237e-05, 
#                        'ch5': 1.166748190444680592e-04,
#                        'ch7': 1.121816242516757850e-04,
#                        'ch9': 1.249543278366777640e-04}
#
#expected_differences = {'Boundary':  8.373150808730501615e-04,
#                        'ch5': 3.425914311580337875e-03,
#                        'ch7': 2.802327594773105189e-03,
#                        'ch9': 3.198733498646373370e-03}


expected_covariances = {'Boundary': 5.299487474489660856e-05,
                        'ch5': 1.169650160375980370e-04,
                        'ch7': 1.123947836450141360e-04,
                        'ch9': 1.248329330436513066e-04}

expected_differences = {'Boundary': 8.269932106586290379e-04,
                        'ch5': 3.411769502897898498e-03, 
                        'ch7': 2.777024544282339583e-03,
                        'ch9': 3.183530359567784788e-03}





#-------------------------
# Read validation dataa
#-------------------------

print 'Reading', project.boundary_filename
fid = NetCDFFile(project.boundary_filename, 'r')
input_time = fid.variables['time'][:]
validation_data['Boundary'] = fid.variables['stage'][:]

reference_time = []
fid = open(project.validation_filename)
lines = fid.readlines()
fid.close()

for i, line in enumerate(lines[1:]):
    if i == len(input_time): break
    
    fields = line.split()

    reference_time.append(float(fields[0]))    # Record reference time
    for j, key in enumerate(gauge_names[1:]):  # Omit boundary gauge
        value = float(fields[1:][j])           # Omit time 
        validation_data[key].append(value/100) # Convert cm2m


# Checks
assert reference_time[0] == 0.0
assert reference_time[-1] == finaltime
assert allclose(reference_time, input_time)

for key in gauge_names:
    validation_data[key] = ensure_numeric(validation_data[key])


#--------------------------------------------------
# Read and interpolate model output
#--------------------------------------------------

import sys
if len(sys.argv) > 1:
    sww_filename = sys.argv[1]
else:    
    sww_filename = project.output_filename
    
#f = file_function('okushiri_new_limiters.sww',   #The best so far
#f = file_function('okushiri_as2005_with_mxspd=0.1.sww',
f = file_function(sww_filename,
                  quantities='stage',
                  interpolation_points=gauge_locations,
                  use_cache=True,
                  verbose=True)


#--------------------------------------------------
# Compare model output to validation data
#--------------------------------------------------

eps = get_machine_precision()
for k, name in enumerate(gauge_names):
    sqsum = 0
    denom = 0
    model = []
    print 
    print 'Validating ' + name
    observed_timeseries = validation_data[name]
    for i, t in enumerate(reference_time):
        model.append(f(t, point_id=k)[0])

    # Covariance measures    
    res = cov(observed_timeseries, model)     
    print 'Covariance = %.18e' %res

    if res < expected_covariances[name]-eps:
        print 'Result is better than expected by: %.18e'\
              %(res-expected_covariances[name])
        print 'Expect = %.18e' %expected_covariances[name]    
    elif res > expected_covariances[name]+eps:
        print 'FAIL: Result is worse than expected by: %.18e'\
                                %(res-expected_covariances[name])
        print 'Expect = %.18e' %expected_covariances[name]
    else:
        pass

    # Difference measures    
    res = sum(abs(observed_timeseries-model))/len(model)     
    print 'Accumulated difference = %.18e' %res

    if res < expected_differences[name]-eps:
        print 'Result is better than expected by: %.18e'\
              %(res-expected_differences[name])
        print 'Expect = %.18e' %expected_differences[name]    
    elif res > expected_differences[name]+eps:
        print 'FAIL: Result is worse than expected by: %.18e'\
                                %(res-expected_differences[name])
        print 'Expect = %.18e' %expected_differences[name]
    else:
        pass




    if plotting is True:
        ion()
        hold(False)
    
        plot(reference_time, validation_data[name], 'r.-',
             reference_time, model, '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')