source: anuga_work/development/stochastic_study/extract_timeseries.py @ 5175

Last change on this file since 5175 was 3514, checked in by duncan, 19 years ago

Hi all,
I'm doing a change in the anuga structure, moving the code to

\anuga_core\source\anuga

After you have done an svn update, the PYTHONPATH has to be changed to;
PYTHONPATH = anuga_core/source/

This is part of changes required to make installation of anuga quicker and reducing the size of our sandpits.

If any imports are broken, try fixing them. With adding anuga. to them for example. If this seems to have really broken things, email/phone me.

Cheers
Duncan

File size: 2.6 KB
Line 
1"""Read in tms file, interpolate at specified locations (ch 5,7,9) and compare
2"""
3
4import sys
5from os import sep
6from caching import cache
7import project
8
9sys.path.append('..'+sep+'..'+sep)
10
11
12gauges = [[0.000, 1.696]]  #Boundary
13gauges += [[4.521, 1.196],  [4.521, 1.696],  [4.521, 2.196]] #Ch 5-7-9
14
15gauge_names = ['Boundary', 'ch5', 'ch7', 'ch9']
16
17
18finaltime = 22.5
19timestep = 0.05
20
21#Read reference data
22
23#Input wave
24filename = 'input_wave.tms'
25print 'Reading', filename
26from Scientific.IO.NetCDF import NetCDFFile
27from Numeric import allclose
28fid = NetCDFFile(filename, 'r')#
29
30input_time = fid.variables['time'][:]
31input_stage = fid.variables['stage'][:]
32
33
34#gauges
35reference_time = []
36ch5 = []
37ch7 = []
38ch9 = []
39filename = 'output_ch5-7-9.txt'
40fid = open(filename)
41lines = fid.readlines()
42fid.close()
43for i, line in enumerate(lines[1:]):
44    if i == len(input_time): break
45
46    fields = line.split()
47
48    reference_time.append(float(fields[0]))
49    ch5.append(float(fields[1])/100)   #cm2m
50    ch7.append(float(fields[2])/100)   #cm2m
51    ch9.append(float(fields[3])/100)   #cm2m
52
53
54from anuga.pyvolution.util import file_function
55from anuga.utilities.numerical_tools import ensure_numeric
56gauge_values = [ensure_numeric(input_stage),
57                ensure_numeric(ch5),
58                ensure_numeric(ch7),
59                ensure_numeric(ch9)] #Reference values
60
61
62
63#Read model output
64#filename = project.basename + '_original.sww'
65filename = project.basename + '.sww'
66
67f = cache(file_function, filename,
68          {'quantities': 'stage',
69           'interpolation_points': gauges,
70           'verbose': True},
71          #evaluate = True,
72          dependencies = [filename],
73          verbose = True)
74
75
76
77
78#Checks
79#print reference_time
80#print input_time
81assert reference_time[0] == 0.0
82assert reference_time[-1] == finaltime
83assert allclose(reference_time, input_time)
84
85
86
87#Validation
88
89
90for k, name in enumerate(gauge_names):
91    sqsum = 0
92    denom = 0
93    model = []
94    print 'Validating ' + name
95    for i, t in enumerate(reference_time):
96        ref = gauge_values[k][i]
97        val = f(t, point_id = k)[0]
98        model.append(val)
99
100        sqsum += (ref - val)**2
101        denom += ref**2
102
103    print sqsum
104    print sqsum/denom
105
106    from pylab import *
107    ion()
108    hold(False)
109    plot(reference_time, gauge_values[k], 'r.-',
110         reference_time, model, 'k-')
111    title('Gauge %s' %name)
112    xlabel('time(s)')
113    ylabel('stage (m)')   
114    legend(('Observed', 'Modelled'), shadow=True, loc='upper left')
115    savefig(name, dpi = 300)
116
117    raw_input('Next')
118show()
119
120
121
122#from pylab import *
123#plot(time, stage)
124#show()
Note: See TracBrowser for help on using the repository browser.