source: inundation/pyvolution/test_util.py @ 1995

Last change on this file since 1995 was 1932, checked in by ole, 19 years ago

First use of Python warning system
File size: 32.7 KB
Line 
1#!/usr/bin/env python
2
3
4import unittest
5from Numeric import zeros, array, allclose
6from math import sqrt, pi
7
8from util import *
9from config import epsilon
10from data_manager import timefile2netcdf
11
12
13def test_function(x, y):
14    return x+y
15
16class Test_Util(unittest.TestCase):
17    def setUp(self):
18        pass
19
20    def tearDown(self):
21        pass
22
23
24    def test_gradient(self):
25        x0 = 0.0; y0 = 0.0; z0 = 0.0
26        x1 = 1.0; y1 = 0.0; z1 = -1.0
27        x2 = 0.0; y2 = 1.0; z2 = 0.0
28
29        zx, zy = gradient(x0, y0, x1, y1, x2, y2, z0, z1, z2)
30
31        assert zx == -1.0
32        assert zy == 0.0
33
34    def test_gradient_more(self):
35        x0 = 2.0/3; y0 = 2.0/3
36        x1=  8.0/3; y1 = 2.0/3
37        x2 = 2.0/3; y2 = 8.0/3
38
39        q0 = 2.0+2.0/3
40        q1 = 8.0+2.0/3
41        q2 = 2.0+8.0/3
42
43        #Gradient of fitted pwl surface
44        a, b = gradient(x0, y0, x1, y1, x2, y2, q0, q1, q2)
45
46        assert abs(a - 3.0) < epsilon
47        assert abs(b - 1.0) < epsilon
48
49
50    def test_gradient2(self):
51        """Test two-point gradient
52        """
53       
54        x0 = 5.0; y0 = 5.0; z0 = 10.0
55        x1 = 8.0; y1 = 2.0; z1 = 1.0
56        x2 = 8.0; y2 = 8.0; z2 = 10.0
57
58        #Reference
59        zx, zy = gradient(x0, y0, x1, y1, x2, y2, z0, z1, z2)
60        a, b = gradient2(x0, y0, x1, y1, z0, z1)
61
62        assert zx == a
63        assert zy == b
64
65        z2_computed = z0 + a*(x2-x0) + b*(y2-y0)
66        assert z2_computed == z2
67       
68    def test_gradient2_more(self):
69        """Test two-point gradient more
70        """
71        x0 = 2.0; y0 = 2.0
72        x1 = 8.0; y1 = 3.0
73        x2 = 1.0; y2 = 8.0
74
75        q0 = 2.0
76        q1 = 8.0
77        q2 = q0
78
79        #Gradient of fitted pwl surface
80        a_ref, b_ref = gradient(x0, y0, x1, y1, x2, y2, q0, q1, q2)
81        a, b = gradient2(x0, y0, x1, y1, q0, q1)       
82
83        assert a == a_ref
84        assert b == b_ref
85
86
87    def test_that_C_extension_compiles(self):
88        FN = 'util_ext.c'
89        try:
90            import util_ext
91        except:
92            from compile import compile
93
94            try:
95                compile(FN)
96            except:
97                raise 'Could not compile %s' %FN
98            else:
99                import util_ext
100
101
102    def test_gradient_C_extension(self):
103        from util_ext import gradient as gradient_c
104
105        x0 = 2.0/3; y0 = 2.0/3
106        x1=  8.0/3; y1 = 2.0/3
107        x2 = 2.0/3; y2 = 8.0/3
108
109        q0 = 2.0+2.0/3
110        q1 = 8.0+2.0/3
111        q2 = 2.0+8.0/3
112
113        #Gradient of fitted pwl surface
114        a, b = gradient_c(x0, y0, x1, y1, x2, y2, q0, q1, q2)
115
116        assert abs(a - 3.0) < epsilon
117        assert abs(b - 1.0) < epsilon
118
119
120    def test_gradient_C_extension3(self):
121        from util_ext import gradient as gradient_c
122
123        from RandomArray import uniform, seed
124        seed(17, 53)
125
126        x0, x1, x2, y0, y1, y2 = uniform(0.0,3.0,6)
127
128        q0 = uniform(0.0, 10.0, 4)
129        q1 = uniform(1.0, 3.0, 4)
130        q2 = uniform(7.0, 20.0, 4)
131
132
133        for i in range(4):
134            #Gradient of fitted pwl surface
135            from util import gradient_python
136            a_ref, b_ref = gradient(x0, y0, x1, y1, x2, y2,
137                                    q0[i], q1[i], q2[i])
138
139            #print a_ref, b_ref
140            a, b = gradient_c(x0, y0, x1, y1, x2, y2,
141                              q0[i], q1[i], q2[i])
142
143            #print a, a_ref, b, b_ref
144            assert abs(a - a_ref) < epsilon
145            assert abs(b - b_ref) < epsilon
146
147
148
149    #Geometric
150    #def test_distance(self):
151    #    from util import distance#
152    #
153    #    self.failUnless( distance([4,2],[7,6]) == 5.0,
154    #                     'Distance is wrong!')
155    #    self.failUnless( allclose(distance([7,6],[9,8]), 2.82842712475),
156    #                    'distance is wrong!')
157    #    self.failUnless( allclose(distance([9,8],[4,2]), 7.81024967591),
158    #                    'distance is wrong!')
159    #
160    #    self.failUnless( distance([9,8],[4,2]) == distance([4,2],[9,8]),
161    #                    'distance is wrong!')
162
163    def test_angle(self):
164        from util import angle
165
166        assert allclose(angle([1.0, 1.0])/pi*180, 45.0)
167
168
169    def test_anglediff(self):
170        from util import anglediff
171
172        assert allclose(anglediff([0.0, 1.], [1.0, 1.0])/pi*180, 45.0)
173
174
175
176    def test_file_function_time1(self):
177        """Test that File function interpolates correctly
178        between given times. No x,y dependency here.
179        """
180
181        #Write file
182        import os, time
183        from config import time_format
184        from math import sin, pi
185
186        #Typical ASCII file
187        finaltime = 1200
188        filename = 'test_file_function'
189        fid = open(filename + '.txt', 'w')
190        start = time.mktime(time.strptime('2000', '%Y'))
191        dt = 60  #One minute intervals
192        t = 0.0
193        while t <= finaltime:
194            t_string = time.strftime(time_format, time.gmtime(t+start))
195            fid.write('%s, %f %f %f\n' %(t_string, 2*t, t**2, sin(t*pi/600)))
196            t += dt
197
198        fid.close()
199
200        #Convert ASCII file to NetCDF (Which is what we really like!)
201        timefile2netcdf(filename)
202
203
204        #Create file function from time series
205        F = file_function(filename + '.tms',
206                          quantities = ['Attribute0',
207                                        'Attribute1',
208                                        'Attribute2'])
209       
210        #Now try interpolation
211        for i in range(20):
212            t = i*10
213            q = F(t)
214
215            #Exact linear intpolation
216            assert allclose(q[0], 2*t)
217            if i%6 == 0:
218                assert allclose(q[1], t**2)
219                assert allclose(q[2], sin(t*pi/600))
220
221        #Check non-exact
222
223        t = 90 #Halfway between 60 and 120
224        q = F(t)
225        assert allclose( (120**2 + 60**2)/2, q[1] )
226        assert allclose( (sin(120*pi/600) + sin(60*pi/600))/2, q[2] )
227
228
229        t = 100 #Two thirds of the way between between 60 and 120
230        q = F(t)
231        assert allclose( 2*120**2/3 + 60**2/3, q[1] )
232        assert allclose( 2*sin(120*pi/600)/3 + sin(60*pi/600)/3, q[2] )
233
234        os.remove(filename + '.txt')
235        os.remove(filename + '.tms')       
236
237
238       
239    def test_spatio_temporal_file_function(self):
240        """Test that spatio temporal file function performs the correct
241        interpolations in both time and space
242        NetCDF version (x,y,t dependency)       
243        """
244        import time
245
246        #Create sww file of simple propagation from left to right
247        #through rectangular domain
248        from shallow_water import Domain, Dirichlet_boundary
249        from mesh_factory import rectangular
250        from Numeric import take, concatenate, reshape
251
252        #Create basic mesh and shallow water domain
253        points, vertices, boundary = rectangular(3, 3)
254        domain1 = Domain(points, vertices, boundary)
255
256        from util import mean
257        domain1.reduction = mean
258        domain1.smooth = True #NOTE: Mimic sww output where each vertex has
259                              # only one value.
260
261        domain1.default_order = 2
262        domain1.store = True
263        domain1.set_datadir('.')
264        domain1.set_name('spatio_temporal_boundary_source_%d' %(id(self)))
265        domain1.quantities_to_be_stored = ['stage', 'xmomentum', 'ymomentum']
266
267        #Bed-slope, friction and IC at vertices (and interpolated elsewhere)
268        domain1.set_quantity('elevation', 0)
269        domain1.set_quantity('friction', 0)
270        domain1.set_quantity('stage', 0)
271
272        # Boundary conditions
273        B0 = Dirichlet_boundary([0,0,0])
274        B6 = Dirichlet_boundary([0.6,0,0])
275        domain1.set_boundary({'left': B6, 'top': B6, 'right': B0, 'bottom': B0})
276        domain1.check_integrity()
277
278        finaltime = 8
279        #Evolution
280        for t in domain1.evolve(yieldstep = 0.1, finaltime = finaltime):
281            pass
282            #domain1.write_time()
283
284
285        #Now read data from sww and check
286        from Scientific.IO.NetCDF import NetCDFFile
287        filename = domain1.get_name() + '.' + domain1.format
288        fid = NetCDFFile(filename)
289
290        x = fid.variables['x'][:]
291        y = fid.variables['y'][:]
292        stage = fid.variables['stage'][:]
293        xmomentum = fid.variables['xmomentum'][:]
294        ymomentum = fid.variables['ymomentum'][:]
295        time = fid.variables['time'][:]
296
297        #Take stage vertex values at last timestep on diagonal
298        #Diagonal is identified by vertices: 0, 5, 10, 15
299
300        timestep = len(time)-1 #Last timestep
301        d_stage = reshape(take(stage[timestep, :], [0,5,10,15]), (4,1))
302        d_uh = reshape(take(xmomentum[timestep, :], [0,5,10,15]), (4,1))
303        d_vh = reshape(take(ymomentum[timestep, :], [0,5,10,15]), (4,1))
304        D = concatenate( (d_stage, d_uh, d_vh), axis=1)
305
306        #Reference interpolated values at midpoints on diagonal at
307        #this timestep are
308        r0 = (D[0] + D[1])/2
309        r1 = (D[1] + D[2])/2
310        r2 = (D[2] + D[3])/2
311
312        #And the midpoints are found now
313        Dx = take(reshape(x, (16,1)), [0,5,10,15])
314        Dy = take(reshape(y, (16,1)), [0,5,10,15])
315
316        diag = concatenate( (Dx, Dy), axis=1)
317        d_midpoints = (diag[1:] + diag[:-1])/2
318
319        #Let us see if the file function can find the correct
320        #values at the midpoints at the last timestep:
321        f = file_function(filename, domain1,
322                          interpolation_points = d_midpoints)
323
324        q = f(timestep/10., point_id=0); assert allclose(r0, q)
325        q = f(timestep/10., point_id=1); assert allclose(r1, q)
326        q = f(timestep/10., point_id=2); assert allclose(r2, q)
327
328
329        ##################
330        #Now do the same for the first timestep
331
332        timestep = 0 #First timestep
333        d_stage = reshape(take(stage[timestep, :], [0,5,10,15]), (4,1))
334        d_uh = reshape(take(xmomentum[timestep, :], [0,5,10,15]), (4,1))
335        d_vh = reshape(take(ymomentum[timestep, :], [0,5,10,15]), (4,1))
336        D = concatenate( (d_stage, d_uh, d_vh), axis=1)
337
338        #Reference interpolated values at midpoints on diagonal at
339        #this timestep are
340        r0 = (D[0] + D[1])/2
341        r1 = (D[1] + D[2])/2
342        r2 = (D[2] + D[3])/2
343
344        #Let us see if the file function can find the correct
345        #values
346        q = f(0, point_id=0); assert allclose(r0, q)
347        q = f(0, point_id=1); assert allclose(r1, q)
348        q = f(0, point_id=2); assert allclose(r2, q)
349
350
351        ##################
352        #Now do it again for a timestep in the middle
353
354        timestep = 33
355        d_stage = reshape(take(stage[timestep, :], [0,5,10,15]), (4,1))
356        d_uh = reshape(take(xmomentum[timestep, :], [0,5,10,15]), (4,1))
357        d_vh = reshape(take(ymomentum[timestep, :], [0,5,10,15]), (4,1))
358        D = concatenate( (d_stage, d_uh, d_vh), axis=1)
359
360        #Reference interpolated values at midpoints on diagonal at
361        #this timestep are
362        r0 = (D[0] + D[1])/2
363        r1 = (D[1] + D[2])/2
364        r2 = (D[2] + D[3])/2
365
366        q = f(timestep/10., point_id=0); assert allclose(r0, q)
367        q = f(timestep/10., point_id=1); assert allclose(r1, q)
368        q = f(timestep/10., point_id=2); assert allclose(r2, q)
369
370
371        ##################
372        #Now check temporal interpolation
373        #Halfway between timestep 15 and 16
374
375        timestep = 15
376        d_stage = reshape(take(stage[timestep, :], [0,5,10,15]), (4,1))
377        d_uh = reshape(take(xmomentum[timestep, :], [0,5,10,15]), (4,1))
378        d_vh = reshape(take(ymomentum[timestep, :], [0,5,10,15]), (4,1))
379        D = concatenate( (d_stage, d_uh, d_vh), axis=1)
380
381        #Reference interpolated values at midpoints on diagonal at
382        #this timestep are
383        r0_0 = (D[0] + D[1])/2
384        r1_0 = (D[1] + D[2])/2
385        r2_0 = (D[2] + D[3])/2
386
387        #
388        timestep = 16
389        d_stage = reshape(take(stage[timestep, :], [0,5,10,15]), (4,1))
390        d_uh = reshape(take(xmomentum[timestep, :], [0,5,10,15]), (4,1))
391        d_vh = reshape(take(ymomentum[timestep, :], [0,5,10,15]), (4,1))
392        D = concatenate( (d_stage, d_uh, d_vh), axis=1)
393
394        #Reference interpolated values at midpoints on diagonal at
395        #this timestep are
396        r0_1 = (D[0] + D[1])/2
397        r1_1 = (D[1] + D[2])/2
398        r2_1 = (D[2] + D[3])/2
399
400        # The reference values are
401        r0 = (r0_0 + r0_1)/2
402        r1 = (r1_0 + r1_1)/2
403        r2 = (r2_0 + r2_1)/2
404
405        q = f((timestep - 0.5)/10., point_id=0); assert allclose(r0, q)
406        q = f((timestep - 0.5)/10., point_id=1); assert allclose(r1, q)
407        q = f((timestep - 0.5)/10., point_id=2); assert allclose(r2, q)
408
409        ##################
410        #Finally check interpolation 2 thirds of the way
411        #between timestep 15 and 16
412
413        # The reference values are
414        r0 = (r0_0 + 2*r0_1)/3
415        r1 = (r1_0 + 2*r1_1)/3
416        r2 = (r2_0 + 2*r2_1)/3
417
418        #And the file function gives
419        q = f((timestep - 1.0/3)/10., point_id=0); assert allclose(r0, q)
420        q = f((timestep - 1.0/3)/10., point_id=1); assert allclose(r1, q)
421        q = f((timestep - 1.0/3)/10., point_id=2); assert allclose(r2, q)
422
423        fid.close()
424        import os
425        os.remove(filename)
426
427
428
429    def test_spatio_temporal_file_function_different_origin(self):
430        """Test that spatio temporal file function performs the correct
431        interpolations in both time and space where space is offset by
432        xllcorner and yllcorner
433        NetCDF version (x,y,t dependency)       
434        """
435        import time
436
437        #Create sww file of simple propagation from left to right
438        #through rectangular domain
439        from shallow_water import Domain, Dirichlet_boundary
440        from mesh_factory import rectangular
441        from Numeric import take, concatenate, reshape
442
443
444        from coordinate_transforms.geo_reference import Geo_reference
445        xllcorner = 2048
446        yllcorner = 11000
447        zone = 2
448
449        #Create basic mesh and shallow water domain
450        points, vertices, boundary = rectangular(3, 3)
451        domain1 = Domain(points, vertices, boundary,
452                         geo_reference = Geo_reference(xllcorner = xllcorner,
453                                                       yllcorner = yllcorner))
454       
455
456        from util import mean
457        domain1.reduction = mean
458        domain1.smooth = True #NOTE: Mimic sww output where each vertex has
459                              # only one value.
460
461        domain1.default_order = 2
462        domain1.store = True
463        domain1.set_datadir('.')
464        domain1.set_name('spatio_temporal_boundary_source_%d' %(id(self)))
465        domain1.quantities_to_be_stored = ['stage', 'xmomentum', 'ymomentum']
466
467        #Bed-slope, friction and IC at vertices (and interpolated elsewhere)
468        domain1.set_quantity('elevation', 0)
469        domain1.set_quantity('friction', 0)
470        domain1.set_quantity('stage', 0)
471
472        # Boundary conditions
473        B0 = Dirichlet_boundary([0,0,0])
474        B6 = Dirichlet_boundary([0.6,0,0])
475        domain1.set_boundary({'left': B6, 'top': B6, 'right': B0, 'bottom': B0})
476        domain1.check_integrity()
477
478        finaltime = 8
479        #Evolution
480        for t in domain1.evolve(yieldstep = 0.1, finaltime = finaltime):
481            pass
482            #domain1.write_time()
483
484
485        #Now read data from sww and check
486        from Scientific.IO.NetCDF import NetCDFFile
487        filename = domain1.get_name() + '.' + domain1.format
488        fid = NetCDFFile(filename)
489
490        x = fid.variables['x'][:]
491        y = fid.variables['y'][:]
492        stage = fid.variables['stage'][:]
493        xmomentum = fid.variables['xmomentum'][:]
494        ymomentum = fid.variables['ymomentum'][:]
495        time = fid.variables['time'][:]
496
497        #Take stage vertex values at last timestep on diagonal
498        #Diagonal is identified by vertices: 0, 5, 10, 15
499
500        timestep = len(time)-1 #Last timestep
501        d_stage = reshape(take(stage[timestep, :], [0,5,10,15]), (4,1))
502        d_uh = reshape(take(xmomentum[timestep, :], [0,5,10,15]), (4,1))
503        d_vh = reshape(take(ymomentum[timestep, :], [0,5,10,15]), (4,1))
504        D = concatenate( (d_stage, d_uh, d_vh), axis=1)
505
506        #Reference interpolated values at midpoints on diagonal at
507        #this timestep are
508        r0 = (D[0] + D[1])/2
509        r1 = (D[1] + D[2])/2
510        r2 = (D[2] + D[3])/2
511
512        #And the midpoints are found now
513        Dx = take(reshape(x, (16,1)), [0,5,10,15])
514        Dy = take(reshape(y, (16,1)), [0,5,10,15])
515
516        diag = concatenate( (Dx, Dy), axis=1)
517        d_midpoints = (diag[1:] + diag[:-1])/2
518
519
520        #Adjust for georef - make interpolation points absolute
521        d_midpoints[:,0] += xllcorner
522        d_midpoints[:,1] += yllcorner               
523
524        #Let us see if the file function can find the correct
525        #values at the midpoints at the last timestep:
526        f = file_function(filename, domain1,
527                          interpolation_points = d_midpoints)
528
529        q = f(timestep/10., point_id=0); assert allclose(r0, q)
530        q = f(timestep/10., point_id=1); assert allclose(r1, q)
531        q = f(timestep/10., point_id=2); assert allclose(r2, q)
532
533
534        ##################
535        #Now do the same for the first timestep
536
537        timestep = 0 #First timestep
538        d_stage = reshape(take(stage[timestep, :], [0,5,10,15]), (4,1))
539        d_uh = reshape(take(xmomentum[timestep, :], [0,5,10,15]), (4,1))
540        d_vh = reshape(take(ymomentum[timestep, :], [0,5,10,15]), (4,1))
541        D = concatenate( (d_stage, d_uh, d_vh), axis=1)
542
543        #Reference interpolated values at midpoints on diagonal at
544        #this timestep are
545        r0 = (D[0] + D[1])/2
546        r1 = (D[1] + D[2])/2
547        r2 = (D[2] + D[3])/2
548
549        #Let us see if the file function can find the correct
550        #values
551        q = f(0, point_id=0); assert allclose(r0, q)
552        q = f(0, point_id=1); assert allclose(r1, q)
553        q = f(0, point_id=2); assert allclose(r2, q)
554
555
556        ##################
557        #Now do it again for a timestep in the middle
558
559        timestep = 33
560        d_stage = reshape(take(stage[timestep, :], [0,5,10,15]), (4,1))
561        d_uh = reshape(take(xmomentum[timestep, :], [0,5,10,15]), (4,1))
562        d_vh = reshape(take(ymomentum[timestep, :], [0,5,10,15]), (4,1))
563        D = concatenate( (d_stage, d_uh, d_vh), axis=1)
564
565        #Reference interpolated values at midpoints on diagonal at
566        #this timestep are
567        r0 = (D[0] + D[1])/2
568        r1 = (D[1] + D[2])/2
569        r2 = (D[2] + D[3])/2
570
571        q = f(timestep/10., point_id=0); assert allclose(r0, q)
572        q = f(timestep/10., point_id=1); assert allclose(r1, q)
573        q = f(timestep/10., point_id=2); assert allclose(r2, q)
574
575
576        ##################
577        #Now check temporal interpolation
578        #Halfway between timestep 15 and 16
579
580        timestep = 15
581        d_stage = reshape(take(stage[timestep, :], [0,5,10,15]), (4,1))
582        d_uh = reshape(take(xmomentum[timestep, :], [0,5,10,15]), (4,1))
583        d_vh = reshape(take(ymomentum[timestep, :], [0,5,10,15]), (4,1))
584        D = concatenate( (d_stage, d_uh, d_vh), axis=1)
585
586        #Reference interpolated values at midpoints on diagonal at
587        #this timestep are
588        r0_0 = (D[0] + D[1])/2
589        r1_0 = (D[1] + D[2])/2
590        r2_0 = (D[2] + D[3])/2
591
592        #
593        timestep = 16
594        d_stage = reshape(take(stage[timestep, :], [0,5,10,15]), (4,1))
595        d_uh = reshape(take(xmomentum[timestep, :], [0,5,10,15]), (4,1))
596        d_vh = reshape(take(ymomentum[timestep, :], [0,5,10,15]), (4,1))
597        D = concatenate( (d_stage, d_uh, d_vh), axis=1)
598
599        #Reference interpolated values at midpoints on diagonal at
600        #this timestep are
601        r0_1 = (D[0] + D[1])/2
602        r1_1 = (D[1] + D[2])/2
603        r2_1 = (D[2] + D[3])/2
604
605        # The reference values are
606        r0 = (r0_0 + r0_1)/2
607        r1 = (r1_0 + r1_1)/2
608        r2 = (r2_0 + r2_1)/2
609
610        q = f((timestep - 0.5)/10., point_id=0); assert allclose(r0, q)
611        q = f((timestep - 0.5)/10., point_id=1); assert allclose(r1, q)
612        q = f((timestep - 0.5)/10., point_id=2); assert allclose(r2, q)
613
614        ##################
615        #Finally check interpolation 2 thirds of the way
616        #between timestep 15 and 16
617
618        # The reference values are
619        r0 = (r0_0 + 2*r0_1)/3
620        r1 = (r1_0 + 2*r1_1)/3
621        r2 = (r2_0 + 2*r2_1)/3
622
623        #And the file function gives
624        q = f((timestep - 1.0/3)/10., point_id=0); assert allclose(r0, q)
625        q = f((timestep - 1.0/3)/10., point_id=1); assert allclose(r1, q)
626        q = f((timestep - 1.0/3)/10., point_id=2); assert allclose(r2, q)
627
628        fid.close()
629        import os
630        os.remove(filename)
631
632       
633
634
635    def test_spatio_temporal_file_function_time(self):
636        """Test that File function interpolates correctly
637        between given times.
638        NetCDF version (x,y,t dependency)
639        """
640
641        #Create NetCDF (sww) file to be read
642        # x: 0, 5, 10, 15
643        # y: -20, -10, 0, 10
644        # t: 0, 60, 120, ...., 1200
645        #
646        # test quantities (arbitrary but non-trivial expressions):
647        #
648        #   stage     = 3*x - y**2 + 2*t
649        #   xmomentum = exp( -((x-7)**2 + (y+5)**2)/20 ) * t**2
650        #   ymomentum = x**2 + y**2 * sin(t*pi/600)
651
652        #NOTE: Nice test that may render some of the others redundant.
653
654        import os, time
655        from config import time_format
656        from Numeric import sin, pi, exp
657        from mesh_factory import rectangular
658        from shallow_water import Domain
659        import data_manager
660
661        finaltime = 1200
662        filename = 'test_file_function'
663
664        #Create a domain to hold test grid
665        #(0:15, -20:10)
666        points, vertices, boundary =\
667                rectangular(4, 4, 15, 30, origin = (0, -20))
668
669
670        #print 'Number of elements', len(vertices)
671        domain = Domain(points, vertices, boundary)
672        domain.smooth = False
673        domain.default_order = 2
674        domain.set_datadir('.')
675        domain.set_name(filename)
676        domain.store = True
677        domain.format = 'sww'   #Native netcdf visualisation format
678
679        #print points
680        start = time.mktime(time.strptime('2000', '%Y'))
681        domain.starttime = start
682
683
684        #Store structure
685        domain.initialise_storage()
686
687        #Compute artificial time steps and store
688        dt = 60  #One minute intervals
689        t = 0.0
690        while t <= finaltime:
691            #Compute quantities
692            f1 = lambda x,y: 3*x - y**2 + 2*t + 4
693            domain.set_quantity('stage', f1)
694
695            f2 = lambda x,y: x+y+t**2
696            domain.set_quantity('xmomentum', f2)
697
698            f3 = lambda x,y: x**2 + y**2 * sin(t*pi/600)
699            domain.set_quantity('ymomentum', f3)
700
701            #Store and advance time
702            domain.time = t
703            domain.store_timestep(domain.conserved_quantities)
704            t += dt
705
706
707        interpolation_points = [[0,-20], [1,0], [0,1], [1.1, 3.14], [10,-12.5]]
708
709
710
711        #Deliberately set domain.starttime to too early
712        domain.starttime = start - 1
713
714        #Create file function
715        F = file_function(filename + '.sww', domain,
716                          quantities = domain.conserved_quantities,
717                          interpolation_points = interpolation_points)
718
719        #Check that FF updates fixes domain starttime
720        assert allclose(domain.starttime, start)
721
722        #Check that domain.starttime isn't updated if later
723        domain.starttime = start + 1
724        F = file_function(filename + '.sww', domain,
725                          quantities = domain.conserved_quantities,
726                          interpolation_points = interpolation_points)
727        assert allclose(domain.starttime, start+1)
728        domain.starttime = start
729
730
731        #Check linear interpolation in time
732        F = file_function(filename + '.sww', domain,
733                          quantities = domain.conserved_quantities,
734                          interpolation_points = interpolation_points)               
735        for id in range(len(interpolation_points)):
736            x = interpolation_points[id][0]
737            y = interpolation_points[id][1]
738
739            for i in range(20):
740                t = i*10
741                k = i%6
742
743                if k == 0:
744                    q0 = F(t, point_id=id)
745                    q1 = F(t+60, point_id=id)
746
747
748                q = F(t, point_id=id)
749                #print i, k, t, q
750                #print ' ', q0
751                #print ' ', q1
752                #print 's', (k*q1 + (6-k)*q0)/6
753                #print
754                assert allclose(q, (k*q1 + (6-k)*q0)/6)
755
756
757        #Another check of linear interpolation in time
758        for id in range(len(interpolation_points)):
759            q60 = F(60, point_id=id)
760            q120 = F(120, point_id=id)
761
762            t = 90 #Halfway between 60 and 120
763            q = F(t, point_id=id)
764            assert allclose( (q120+q60)/2, q )
765
766            t = 100 #Two thirds of the way between between 60 and 120
767            q = F(t, point_id=id)
768            assert allclose(q60/3 + 2*q120/3, q)
769
770
771
772        #Check that domain.starttime isn't updated if later than file starttime but earlier
773        #than file end time
774        delta = 23
775        domain.starttime = start + delta
776        F = file_function(filename + '.sww', domain,
777                          quantities = domain.conserved_quantities,
778                          interpolation_points = interpolation_points)
779        assert allclose(domain.starttime, start+delta)
780
781
782
783
784        #Now try interpolation with delta offset
785        for id in range(len(interpolation_points)):           
786            x = interpolation_points[id][0]
787            y = interpolation_points[id][1]
788
789            for i in range(20):
790                t = i*10
791                k = i%6
792
793                if k == 0:
794                    q0 = F(t-delta, point_id=id)
795                    q1 = F(t+60-delta, point_id=id)
796
797                q = F(t-delta, point_id=id)
798                assert allclose(q, (k*q1 + (6-k)*q0)/6)
799
800
801        os.remove(filename + '.sww')
802
803
804
805    def test_file_function_time_with_domain(self):
806        """Test that File function interpolates correctly
807        between given times. No x,y dependency here.
808        Use domain with starttime
809        """
810
811        #Write file
812        import os, time, calendar
813        from config import time_format
814        from math import sin, pi
815        from domain import Domain
816
817        finaltime = 1200
818        filename = 'test_file_function'
819        fid = open(filename + '.txt', 'w')
820        start = time.mktime(time.strptime('2000', '%Y'))
821        dt = 60  #One minute intervals
822        t = 0.0
823        while t <= finaltime:
824            t_string = time.strftime(time_format, time.gmtime(t+start))
825            fid.write('%s, %f %f %f\n' %(t_string, 2*t, t**2, sin(t*pi/600)))
826            t += dt
827
828        fid.close()
829
830
831        #Convert ASCII file to NetCDF (Which is what we really like!)
832        timefile2netcdf(filename)
833
834
835
836        a = [0.0, 0.0]
837        b = [4.0, 0.0]
838        c = [0.0, 3.0]
839
840        points = [a, b, c]
841        vertices = [[0,1,2]]
842        domain = Domain(points, vertices)
843
844        #Check that domain.starttime is updated if non-existing
845        F = file_function(filename + '.tms', domain)
846       
847        assert allclose(domain.starttime, start)
848
849        #Check that domain.starttime is updated if too early
850        domain.starttime = start - 1
851        F = file_function(filename + '.tms', domain)
852        assert allclose(domain.starttime, start)
853
854        #Check that domain.starttime isn't updated if later
855        domain.starttime = start + 1
856        F = file_function(filename + '.tms', domain)       
857        assert allclose(domain.starttime, start+1)
858
859        domain.starttime = start
860        F = file_function(filename + '.tms', domain,
861                          quantities = ['Attribute0', 'Attribute1', 'Attribute2'])
862       
863
864        #print F.T
865        #print F.precomputed_values
866        #print 'F(60)', F(60)
867       
868        #Now try interpolation
869        for i in range(20):
870            t = i*10
871            q = F(t)
872
873            #Exact linear intpolation
874            assert allclose(q[0], 2*t)
875            if i%6 == 0:
876                assert allclose(q[1], t**2)
877                assert allclose(q[2], sin(t*pi/600))
878
879        #Check non-exact
880
881        t = 90 #Halfway between 60 and 120
882        q = F(t)
883        assert allclose( (120**2 + 60**2)/2, q[1] )
884        assert allclose( (sin(120*pi/600) + sin(60*pi/600))/2, q[2] )
885
886
887        t = 100 #Two thirds of the way between between 60 and 120
888        q = F(t)
889        assert allclose( 2*120**2/3 + 60**2/3, q[1] )
890        assert allclose( 2*sin(120*pi/600)/3 + sin(60*pi/600)/3, q[2] )
891
892        os.remove(filename + '.tms')
893        os.remove(filename + '.txt')       
894
895    def test_file_function_time_with_domain_different_start(self):
896        """Test that File function interpolates correctly
897        between given times. No x,y dependency here.
898        Use domain with a starttime later than that of file
899
900        ASCII version
901        """
902
903        #Write file
904        import os, time, calendar
905        from config import time_format
906        from math import sin, pi
907        from domain import Domain
908
909        finaltime = 1200
910        filename = 'test_file_function'
911        fid = open(filename + '.txt', 'w')
912        start = time.mktime(time.strptime('2000', '%Y'))
913        dt = 60  #One minute intervals
914        t = 0.0
915        while t <= finaltime:
916            t_string = time.strftime(time_format, time.gmtime(t+start))
917            fid.write('%s, %f %f %f\n' %(t_string, 2*t, t**2, sin(t*pi/600)))
918            t += dt
919
920        fid.close()
921
922        #Convert ASCII file to NetCDF (Which is what we really like!)
923        timefile2netcdf(filename)       
924
925        a = [0.0, 0.0]
926        b = [4.0, 0.0]
927        c = [0.0, 3.0]
928
929        points = [a, b, c]
930        vertices = [[0,1,2]]
931        domain = Domain(points, vertices)
932
933        #Check that domain.starttime isn't updated if later than file starttime but earlier
934        #than file end time
935        delta = 23
936        domain.starttime = start + delta
937        F = file_function(filename + '.tms', domain,
938                          quantities = ['Attribute0', 'Attribute1', 'Attribute2'])       
939        assert allclose(domain.starttime, start+delta)
940
941
942
943
944        #Now try interpolation with delta offset
945        for i in range(20):
946            t = i*10
947            q = F(t-delta)
948
949            #Exact linear intpolation
950            assert allclose(q[0], 2*t)
951            if i%6 == 0:
952                assert allclose(q[1], t**2)
953                assert allclose(q[2], sin(t*pi/600))
954
955        #Check non-exact
956
957        t = 90 #Halfway between 60 and 120
958        q = F(t-delta)
959        assert allclose( (120**2 + 60**2)/2, q[1] )
960        assert allclose( (sin(120*pi/600) + sin(60*pi/600))/2, q[2] )
961
962
963        t = 100 #Two thirds of the way between between 60 and 120
964        q = F(t-delta)
965        assert allclose( 2*120**2/3 + 60**2/3, q[1] )
966        assert allclose( 2*sin(120*pi/600)/3 + sin(60*pi/600)/3, q[2] )
967
968
969        os.remove(filename + '.tms')
970        os.remove(filename + '.txt')               
971
972
973
974    def test_apply_expression_to_dictionary(self):
975
976        foo = array([[1,2,3],
977                     [4,5,6]])
978
979        bar = array([[-1,0,5],
980                     [6,1,1]])                 
981
982        D = {'X': foo, 'Y': bar}
983
984        Z = apply_expression_to_dictionary('X+Y', D)       
985        assert allclose(Z, foo+bar)
986
987        Z = apply_expression_to_dictionary('X*Y', D)       
988        assert allclose(Z, foo*bar)       
989
990        Z = apply_expression_to_dictionary('4*X+Y', D)       
991        assert allclose(Z, 4*foo+bar)       
992
993        #Check exceptions
994        try:
995            #Wrong name
996            Z = apply_expression_to_dictionary('4*X+A', D)       
997        except NameError:
998            pass
999        else:
1000            msg = 'Should have raised a NameError Exception'
1001            raise msg
1002
1003
1004        try:
1005            #Wrong order
1006            Z = apply_expression_to_dictionary(D, '4*X+A')       
1007        except AssertionError:
1008            pass
1009        else:
1010            msg = 'Should have raised a AssertionError Exception'
1011            raise msg       
1012       
1013
1014    def test_multiple_replace(self):
1015        """Hard test that checks a true word-by-word simultaneous replace
1016        """
1017       
1018        D = {'x': 'xi', 'y': 'eta', 'xi':'lam'}
1019        exp = '3*x+y + xi'
1020       
1021        new = multiple_replace(exp, D)
1022       
1023        assert new == '3*xi+eta + lam'
1024                         
1025
1026
1027    def test_point_on_line_obsolete(self):
1028        """Test that obsolete call issues appropriate warning"""
1029
1030        #Turn warning into an exception
1031        import warnings
1032        warnings.filterwarnings('error')
1033
1034        try:
1035            assert point_on_line( 0, 0.5, 0,1, 0,0 )
1036        except DeprecationWarning:
1037            pass
1038        else:
1039            msg = 'point_on_line should have issued a DeprecationWarning'
1040            raise Exception(msg)   
1041
1042        warnings.resetwarnings()                         
1043                         
1044
1045#-------------------------------------------------------------
1046if __name__ == "__main__":
1047    suite = unittest.makeSuite(Test_Util,'test')
1048    #suite = unittest.makeSuite(Test_Util,'test_file_function_time')
1049    runner = unittest.TextTestRunner()
1050    runner.run(suite)
1051
1052
1053
1054
Note: See TracBrowser for help on using the repository browser.