source: trunk/anuga_core/source/anuga/shallow_water/test_forcing.py @ 8116

Last change on this file since 8116 was 8116, checked in by jakeman, 13 years ago

jakeman: updated forcing.py to include wind and pressure stress terms. Updated test_forcing.py with associated unit tests
File size: 64.4 KB
Line 
1"""  Test environmental forcing - rain, wind, etc.
2"""
3
4import unittest, os
5import anuga
6from anuga.shallow_water.shallow_water_domain import Domain
7from boundaries import Reflective_boundary
8from anuga.coordinate_transforms.geo_reference import Geo_reference
9from anuga.file_conversion.file_conversion import timefile2netcdf
10from forcing import *
11from mesh_factory import rectangular
12from file_conversion.sts2sww_mesh import sts2sww_mesh
13from anuga.abstract_2d_finite_volumes.util import file_function
14
15import numpy as num
16import warnings
17
18
19def scalar_func_list(t, x, y):
20    """Function that returns a scalar.
21
22    Used to test error message when numeric array is expected
23    """
24
25    return [17.7]
26
27
28def speed(t, x, y):
29    """
30    Variable windfield implemented using functions   
31    Large speeds halfway between center and edges
32
33    Low speeds at center and edges
34    """
35
36    from math import exp, cos, pi
37
38    x = num.array(x)
39    y = num.array(y)
40
41    N = len(x)
42    s = 0*x  #New array
43
44    for k in range(N):
45        r = num.sqrt(x[k]**2 + y[k]**2)
46        factor = exp(-(r-0.15)**2)
47        s[k] = 4000 * factor * (cos(t*2*pi/150) + 2)
48
49    return s
50
51
52def angle(t, x, y):
53    """Rotating field
54    """
55    from math import atan, pi
56
57    x = num.array(x)
58    y = num.array(y)
59
60    N = len(x)
61    a = 0 * x    # New array
62
63    for k in range(N):
64        r = num.sqrt(x[k]**2 + y[k]**2)
65
66        angle = atan(y[k]/x[k])
67
68        if x[k] < 0:
69            angle += pi
70
71        # Take normal direction
72        angle -= pi/2
73
74        # Ensure positive radians
75        if angle < 0:
76            angle += 2*pi
77
78        a[k] = angle/pi*180
79
80    return a
81
82def time_varying_speed(t, x, y):
83    """
84    Variable speed windfield
85    """
86
87    from math import exp, cos, pi
88
89    x = num.array(x,num.float)
90    y = num.array(y,num.float)
91
92    N = len(x)
93    s = 0*x  #New array
94
95    #dx=x[-1]-x[0]; dy = y[-1]-y[0]
96    S=100.
97    for k in range(N):
98        s[k]=S*(1.+t/100.)
99    return s
100
101
102def time_varying_angle(t, x, y):
103    """Rotating field
104    """
105    from math import atan, pi
106
107    x = num.array(x,num.float)
108    y = num.array(y,num.float)
109
110    N = len(x)
111    a = 0 * x    # New array
112
113    phi=135.
114    for k in range(N):
115        a[k]=phi*(1.+t/100.)
116
117    return a
118
119
120def time_varying_pressure(t, x, y):
121    """Rotating field
122    """
123    from math import atan, pi
124
125    x = num.array(x,num.float)
126    y = num.array(y,num.float)
127
128    N = len(x)
129    p = 0 * x    # New array
130
131    p0=1000.
132    for k in range(N):
133        p[k]=p0*(1.-t/100.)
134
135    return p
136
137def spatial_linear_varying_speed(t, x, y):
138    """
139    Variable speed windfield
140    """
141
142    from math import exp, cos, pi
143
144    x = num.array(x)
145    y = num.array(y)
146
147    N = len(x)
148    s = 0*x  #New array
149
150    #dx=x[-1]-x[0]; dy = y[-1]-y[0]
151    s0=250.
152    ymin=num.min(y)
153    xmin=num.min(x)
154    a=0.000025; b=0.0000125;
155    for k in range(N):
156        s[k]=s0*(1+t/100.)+a*x[k]+b*y[k]
157    return s
158
159
160def spatial_linear_varying_angle(t, x, y):
161    """Rotating field
162    """
163    from math import atan, pi
164
165    x = num.array(x)
166    y = num.array(y)
167
168    N = len(x)
169    a = 0 * x    # New array
170
171    phi=135.
172    b1=0.000025; b2=0.00001125;
173    for k in range(N):
174        a[k]=phi*(1+t/100.)+b1*x[k]+b2*y[k]
175    return a
176
177def spatial_linear_varying_pressure(t, x, y):
178    p0=1000;
179    a=0.000025; b=0.0000125;
180
181    x = num.array(x)
182    y = num.array(y)
183
184    N = len(x)
185    p = 0 * x    # New array
186
187    for k in range(N):
188        p[k]=p0*(1.-t/100.)+a*x[k]+b*y[k]
189    return p
190
191
192def grid_1d(x0,dx,nx):
193    x = num.empty(nx,dtype=num.float)
194    for i in range(nx):
195        x[i]=x0+float(i)*dx
196    return x
197   
198
199def ndgrid(x,y):
200    nx = len(x)
201    ny = len(y)
202    X = num.empty(nx*ny,dtype=num.float)
203    Y = num.empty(nx*ny,dtype=num.float)
204    k=0
205    for i in range(nx):
206        for j in range(ny):
207            X[k]=x[i]
208            Y[k]=y[j]
209            k+=1
210    return X,Y
211
212class Test_Forcing(unittest.TestCase):
213    def setUp(self):
214        pass
215
216    def tearDown(self):
217        pass
218       
219    def write_wind_pressure_field_sts(self,
220                                      field_sts_filename,
221                                      nrows=10,
222                                      ncols=10,
223                                      cellsize=25,
224                                      origin=(0.0,0.0),
225                                      refzone=50,
226                                      timestep=1,
227                                      number_of_timesteps=10,
228                                      angle=135.0,
229                                      speed=100.0,
230                                      pressure=1000.0):
231
232        xllcorner=origin[0]
233        yllcorner=origin[1]
234        starttime = 0; endtime = number_of_timesteps*timestep;
235        no_data = -9999
236
237        time = num.arange(starttime, endtime, timestep, dtype='i')
238
239        x = grid_1d(xllcorner,cellsize,ncols)
240        y = grid_1d(yllcorner,cellsize,nrows)
241        [X,Y] = ndgrid(x,y)
242        number_of_points = nrows*ncols
243
244        wind_speed = num.empty((number_of_timesteps,nrows*ncols),dtype=num.float)
245        wind_angle = num.empty((number_of_timesteps,nrows*ncols),dtype=num.float)
246        barometric_pressure = num.empty((number_of_timesteps,nrows*ncols),
247                                        dtype=num.float)
248
249        if ( callable(speed) and callable(angle) and callable(pressure) ):
250            x = num.ones(3, num.float)
251            y = num.ones(3, num.float)
252            try:
253                s = speed(1.0, x=x, y=y)
254                a = angle(1.0, x=x, y=y)
255                p = pressure(1.0, x=x, y=y)
256                use_function=True
257            except Exception, e:
258                msg = 'Function could not be executed.\n'
259                raise Exception, msg
260        else:
261            try :
262                speed=float(speed)
263                angle=float(angle)
264                pressure=float(pressure)
265                use_function=False
266            except:
267                msg = ('Force fields must be a scalar value coercible to float.')
268                raise Exception, msg
269
270        for i,t in enumerate(time):
271            if ( use_function ):
272                wind_speed[i,:] = speed(t,X,Y)
273                wind_angle[i,:] = angle(t,X,Y)
274                barometric_pressure[i,:] = pressure(t,X,Y)
275            else:
276                wind_speed[i,:] = speed
277                wind_angle[i,:] = angle
278                barometric_pressure[i,:] = pressure
279
280        # "Creating the field STS NetCDF file"
281
282        fid = NetCDFFile(field_sts_filename+'.sts', 'w')
283        fid.institution = 'Geoscience Australia'
284        fid.description = "description"
285        fid.starttime = 0.0
286        fid.ncols = ncols
287        fid.nrows = nrows
288        fid.cellsize = cellsize
289        fid.no_data = no_data
290        fid.createDimension('number_of_points', number_of_points)
291        fid.createDimension('number_of_timesteps', number_of_timesteps)
292        fid.createDimension('numbers_in_range', 2)
293
294        fid.createVariable('x', 'd', ('number_of_points',))
295        fid.createVariable('y', 'd', ('number_of_points',))
296        fid.createVariable('time', 'i', ('number_of_timesteps',))
297        fid.createVariable('wind_speed', 'd', ('number_of_timesteps', 
298                                               'number_of_points'))
299        fid.createVariable('wind_speed_range', 'd', ('numbers_in_range', ))
300        fid.createVariable('wind_angle', 'd', ('number_of_timesteps', 
301                                               'number_of_points'))
302        fid.createVariable('wind_angle_range', 'd', ('numbers_in_range',))
303        fid.createVariable('barometric_pressure', 'd', ('number_of_timesteps', 
304                                             'number_of_points'))
305        fid.createVariable('barometric_pressure_range', 'd', ('numbers_in_range',))
306
307
308        fid.variables['wind_speed_range'][:] = num.array([1e+036, -1e+036])
309        fid.variables['wind_angle_range'][:] = num.array([1e+036, -1e+036])
310        fid.variables['barometric_pressure_range'][:] = num.array([1e+036, -1e+036])
311        fid.variables['time'][:] = time
312
313        ws = fid.variables['wind_speed']
314        wa = fid.variables['wind_angle']
315        pr = fid.variables['barometric_pressure']
316
317        for i in xrange(number_of_timesteps):
318            ws[i] = wind_speed[i,:]
319            wa[i] = wind_angle[i,:]
320            pr[i] = barometric_pressure[i,:]
321
322        origin = anuga.coordinate_transforms.geo_reference.Geo_reference(refzone,
323                                                                         xllcorner,
324                                                                         yllcorner)
325        geo_ref = anuga.coordinate_transforms.geo_reference.write_NetCDF_georeference(origin, fid)
326
327        fid.variables['x'][:]=X-geo_ref.get_xllcorner()
328        fid.variables['y'][:]=Y-geo_ref.get_yllcorner()
329
330
331        fid.close()
332
333    def test_constant_wind_stress(self):
334        from anuga.config import rho_a, rho_w, eta_w
335        from math import pi, cos, sin
336
337        a = [0.0, 0.0]
338        b = [0.0, 2.0]
339        c = [2.0, 0.0]
340        d = [0.0, 4.0]
341        e = [2.0, 2.0]
342        f = [4.0, 0.0]
343
344        points = [a, b, c, d, e, f]
345        #             bac,     bce,     ecf,     dbe
346        vertices = [[1,0,2], [1,2,4], [4,2,5], [3,1,4]]
347
348        domain = Domain(points, vertices)
349
350        #Flat surface with 1m of water
351        domain.set_quantity('elevation', 0)
352        domain.set_quantity('stage', 1.0)
353        domain.set_quantity('friction', 0)
354
355        Br = Reflective_boundary(domain)
356        domain.set_boundary({'exterior': Br})
357
358        #Setup only one forcing term, constant wind stress
359        s = 100
360        phi = 135
361        domain.forcing_terms = []
362        domain.forcing_terms.append(Wind_stress(s, phi))
363
364        domain.compute_forcing_terms()
365
366        const = eta_w*rho_a / rho_w
367
368        #Convert to radians
369        phi = phi*pi / 180
370
371        #Compute velocity vector (u, v)
372        u = s*cos(phi)
373        v = s*sin(phi)
374
375        #Compute wind stress
376        S = const * num.sqrt(u**2 + v**2)
377
378        assert num.allclose(domain.quantities['stage'].explicit_update, 0)
379        assert num.allclose(domain.quantities['xmomentum'].explicit_update, S*u)
380        assert num.allclose(domain.quantities['ymomentum'].explicit_update, S*v)
381
382    def test_variable_wind_stress(self):
383        from anuga.config import rho_a, rho_w, eta_w
384        from math import pi, cos, sin
385
386        a = [0.0, 0.0]
387        b = [0.0, 2.0]
388        c = [2.0, 0.0]
389        d = [0.0, 4.0]
390        e = [2.0, 2.0]
391        f = [4.0, 0.0]
392
393        points = [a, b, c, d, e, f]
394        #             bac,     bce,     ecf,     dbe
395        vertices = [[1,0,2], [1,2,4], [4,2,5], [3,1,4]]
396
397        domain = Domain(points, vertices)
398
399        #Flat surface with 1m of water
400        domain.set_quantity('elevation', 0)
401        domain.set_quantity('stage', 1.0)
402        domain.set_quantity('friction', 0)
403
404        Br = Reflective_boundary(domain)
405        domain.set_boundary({'exterior': Br})
406
407        domain.time = 5.54    # Take a random time (not zero)
408
409        #Setup only one forcing term, constant wind stress
410        s = 100
411        phi = 135
412        domain.forcing_terms = []
413        domain.forcing_terms.append(Wind_stress(s=speed, phi=angle))
414
415        domain.compute_forcing_terms()
416
417        #Compute reference solution
418        const = eta_w*rho_a / rho_w
419
420        N = len(domain)    # number_of_triangles
421
422        xc = domain.get_centroid_coordinates()
423        t = domain.time
424
425        x = xc[:,0]
426        y = xc[:,1]
427        s_vec = speed(t,x,y)
428        phi_vec = angle(t,x,y)
429
430        for k in range(N):
431            # Convert to radians
432            phi = phi_vec[k]*pi / 180
433            s = s_vec[k]
434
435            # Compute velocity vector (u, v)
436            u = s*cos(phi)
437            v = s*sin(phi)
438
439            # Compute wind stress
440            S = const * num.sqrt(u**2 + v**2)
441
442            assert num.allclose(domain.quantities['stage'].explicit_update[k],
443                                0)
444            assert num.allclose(domain.quantities['xmomentum'].\
445                                     explicit_update[k],
446                                S*u)
447            assert num.allclose(domain.quantities['ymomentum'].\
448                                     explicit_update[k],
449                                S*v)
450
451    def test_windfield_from_file(self):
452        import time
453        from anuga.config import rho_a, rho_w, eta_w
454        from math import pi, cos, sin
455        from anuga.config import time_format
456        from anuga.abstract_2d_finite_volumes.util import file_function
457
458        a = [0.0, 0.0]
459        b = [0.0, 2.0]
460        c = [2.0, 0.0]
461        d = [0.0, 4.0]
462        e = [2.0, 2.0]
463        f = [4.0, 0.0]
464
465        points = [a, b, c, d, e, f]
466        #             bac,     bce,     ecf,     dbe
467        vertices = [[1,0,2], [1,2,4], [4,2,5], [3,1,4]]
468
469        domain = Domain(points, vertices)
470
471        # Flat surface with 1m of water
472        domain.set_quantity('elevation', 0)
473        domain.set_quantity('stage', 1.0)
474        domain.set_quantity('friction', 0)
475
476        Br = Reflective_boundary(domain)
477        domain.set_boundary({'exterior': Br})
478
479        domain.time = 7    # Take a time that is represented in file (not zero)
480
481        # Write wind stress file (ensure that domain.time is covered)
482        # Take x=1 and y=0
483        filename = 'test_windstress_from_file'
484        start = time.mktime(time.strptime('2000', '%Y'))
485        fid = open(filename + '.txt', 'w')
486        dt = 1    # One second interval
487        t = 0.0
488        while t <= 10.0:
489            t_string = time.strftime(time_format, time.gmtime(t+start))
490
491            fid.write('%s, %f %f\n' %
492                      (t_string, speed(t,[1],[0])[0], angle(t,[1],[0])[0]))
493            t += dt
494
495        fid.close()
496
497        timefile2netcdf(filename + '.txt')
498        os.remove(filename + '.txt')
499
500        # Setup wind stress
501        F = file_function(filename + '.tms',
502                          quantities=['Attribute0', 'Attribute1'])
503        os.remove(filename + '.tms')
504
505        W = Wind_stress(F)
506
507        domain.forcing_terms = []
508        domain.forcing_terms.append(W)
509
510        domain.compute_forcing_terms()
511
512        # Compute reference solution
513        const = eta_w*rho_a / rho_w
514
515        N = len(domain)    # number_of_triangles
516
517        t = domain.time
518
519        s = speed(t, [1], [0])[0]
520        phi = angle(t, [1], [0])[0]
521
522        # Convert to radians
523        phi = phi*pi / 180
524
525        # Compute velocity vector (u, v)
526        u = s*cos(phi)
527        v = s*sin(phi)
528
529        # Compute wind stress
530        S = const * num.sqrt(u**2 + v**2)
531
532        for k in range(N):
533            assert num.allclose(domain.quantities['stage'].explicit_update[k],
534                                0)
535            assert num.allclose(domain.quantities['xmomentum'].\
536                                    explicit_update[k],
537                                S*u)
538            assert num.allclose(domain.quantities['ymomentum'].\
539                                    explicit_update[k],
540                                S*v)
541
542    def test_windfield_from_file_seconds(self):
543        import time
544        from anuga.config import rho_a, rho_w, eta_w
545        from math import pi, cos, sin
546        from anuga.config import time_format
547        from anuga.abstract_2d_finite_volumes.util import file_function
548
549        a = [0.0, 0.0]
550        b = [0.0, 2.0]
551        c = [2.0, 0.0]
552        d = [0.0, 4.0]
553        e = [2.0, 2.0]
554        f = [4.0, 0.0]
555
556        points = [a, b, c, d, e, f]
557        #             bac,     bce,     ecf,     dbe
558        vertices = [[1,0,2], [1,2,4], [4,2,5], [3,1,4]]
559
560        domain = Domain(points, vertices)
561
562        # Flat surface with 1m of water
563        domain.set_quantity('elevation', 0)
564        domain.set_quantity('stage', 1.0)
565        domain.set_quantity('friction', 0)
566
567        Br = Reflective_boundary(domain)
568        domain.set_boundary({'exterior': Br})
569
570        domain.time = 7    # Take a time that is represented in file (not zero)
571
572        # Write wind stress file (ensure that domain.time is covered)
573        # Take x=1 and y=0
574        filename = 'test_windstress_from_file'
575        start = time.mktime(time.strptime('2000', '%Y'))
576        fid = open(filename + '.txt', 'w')
577        dt = 0.5    # Half second interval
578        t = 0.0
579        while t <= 10.0:
580            fid.write('%s, %f %f\n'
581                      % (str(t), speed(t, [1], [0])[0], angle(t, [1], [0])[0]))
582            t += dt
583
584        fid.close()
585
586        timefile2netcdf(filename + '.txt', time_as_seconds=True)
587        os.remove(filename + '.txt')
588
589        # Setup wind stress
590        F = file_function(filename + '.tms',
591                          quantities=['Attribute0', 'Attribute1'])
592        os.remove(filename + '.tms')
593
594        W = Wind_stress(F)
595
596        domain.forcing_terms = []
597        domain.forcing_terms.append(W)
598
599        domain.compute_forcing_terms()
600
601        # Compute reference solution
602        const = eta_w*rho_a / rho_w
603
604        N = len(domain)    # number_of_triangles
605
606        t = domain.time
607
608        s = speed(t, [1], [0])[0]
609        phi = angle(t, [1], [0])[0]
610
611        # Convert to radians
612        phi = phi*pi / 180
613
614        # Compute velocity vector (u, v)
615        u = s*cos(phi)
616        v = s*sin(phi)
617
618        # Compute wind stress
619        S = const * num.sqrt(u**2 + v**2)
620
621        for k in range(N):
622            assert num.allclose(domain.quantities['stage'].explicit_update[k],
623                                0)
624            assert num.allclose(domain.quantities['xmomentum'].\
625                                    explicit_update[k],
626                                S*u)
627            assert num.allclose(domain.quantities['ymomentum'].\
628                                    explicit_update[k],
629                                S*v)
630
631    def test_wind_stress_error_condition(self):
632        """Test that windstress reacts properly when forcing functions
633        are wrong - e.g. returns a scalar
634        """
635
636        from math import pi, cos, sin
637        from anuga.config import rho_a, rho_w, eta_w
638
639        a = [0.0, 0.0]
640        b = [0.0, 2.0]
641        c = [2.0, 0.0]
642        d = [0.0, 4.0]
643        e = [2.0, 2.0]
644        f = [4.0, 0.0]
645
646        points = [a, b, c, d, e, f]
647        #             bac,     bce,     ecf,     dbe
648        vertices = [[1,0,2], [1,2,4], [4,2,5], [3,1,4]]
649
650        domain = Domain(points, vertices)
651
652        # Flat surface with 1m of water
653        domain.set_quantity('elevation', 0)
654        domain.set_quantity('stage', 1.0)
655        domain.set_quantity('friction', 0)
656
657        Br = Reflective_boundary(domain)
658        domain.set_boundary({'exterior': Br})
659
660        domain.time = 5.54    # Take a random time (not zero)
661
662        # Setup only one forcing term, bad func
663        domain.forcing_terms = []
664
665        try:
666            domain.forcing_terms.append(Wind_stress(s=scalar_func_list,
667                                                    phi=angle))
668        except AssertionError:
669            pass
670        else:
671            msg = 'Should have raised exception'
672            raise Exception, msg
673
674        try:
675            domain.forcing_terms.append(Wind_stress(s=speed, phi=scalar_func))
676        except Exception:
677            pass
678        else:
679            msg = 'Should have raised exception'
680            raise Exception, msg
681
682        try:
683            domain.forcing_terms.append(Wind_stress(s=speed, phi='xx'))
684        except:
685            pass
686        else:
687            msg = 'Should have raised exception'
688            raise Exception, msg
689
690    def test_rainfall(self):
691        from math import pi, cos, sin
692
693        a = [0.0, 0.0]
694        b = [0.0, 2.0]
695        c = [2.0, 0.0]
696        d = [0.0, 4.0]
697        e = [2.0, 2.0]
698        f = [4.0, 0.0]
699
700        points = [a, b, c, d, e, f]
701        #             bac,     bce,     ecf,     dbe
702        vertices = [[1,0,2], [1,2,4], [4,2,5], [3,1,4]]
703
704        domain = Domain(points, vertices)
705
706        # Flat surface with 1m of water
707        domain.set_quantity('elevation', 0)
708        domain.set_quantity('stage', 1.0)
709        domain.set_quantity('friction', 0)
710
711        Br = Reflective_boundary(domain)
712        domain.set_boundary({'exterior': Br})
713
714        # Setup only one forcing term, constant rainfall
715        domain.forcing_terms = []
716        domain.forcing_terms.append(Rainfall(domain, rate=2.0))
717
718        domain.compute_forcing_terms()
719        assert num.allclose(domain.quantities['stage'].explicit_update,
720                            2.0/1000)
721
722    def test_rainfall_restricted_by_polygon(self):
723        from math import pi, cos, sin
724
725        a = [0.0, 0.0]
726        b = [0.0, 2.0]
727        c = [2.0, 0.0]
728        d = [0.0, 4.0]
729        e = [2.0, 2.0]
730        f = [4.0, 0.0]
731
732        points = [a, b, c, d, e, f]
733        #             bac,     bce,     ecf,     dbe
734        vertices = [[1,0,2], [1,2,4], [4,2,5], [3,1,4]]
735
736        domain = Domain(points, vertices)
737
738        # Flat surface with 1m of water
739        domain.set_quantity('elevation', 0)
740        domain.set_quantity('stage', 1.0)
741        domain.set_quantity('friction', 0)
742
743        Br = Reflective_boundary(domain)
744        domain.set_boundary({'exterior': Br})
745
746        # Setup only one forcing term, constant rainfall
747        # restricted to a polygon enclosing triangle #1 (bce)
748        domain.forcing_terms = []
749        R = Rainfall(domain, rate=2.0, polygon=[[1,1], [2,1], [2,2], [1,2]])
750
751        assert num.allclose(R.exchange_area, 2)
752
753        domain.forcing_terms.append(R)
754
755        domain.compute_forcing_terms()
756
757        assert num.allclose(domain.quantities['stage'].explicit_update[1],
758                            2.0/1000)
759        assert num.allclose(domain.quantities['stage'].explicit_update[0], 0)
760        assert num.allclose(domain.quantities['stage'].explicit_update[2:], 0)
761
762    def test_time_dependent_rainfall_restricted_by_polygon(self):
763        a = [0.0, 0.0]
764        b = [0.0, 2.0]
765        c = [2.0, 0.0]
766        d = [0.0, 4.0]
767        e = [2.0, 2.0]
768        f = [4.0, 0.0]
769
770        points = [a, b, c, d, e, f]
771        #             bac,     bce,     ecf,     dbe
772        vertices = [[1,0,2], [1,2,4], [4,2,5], [3,1,4]]
773
774        domain = Domain(points, vertices)
775
776        # Flat surface with 1m of water
777        domain.set_quantity('elevation', 0)
778        domain.set_quantity('stage', 1.0)
779        domain.set_quantity('friction', 0)
780
781        Br = Reflective_boundary(domain)
782        domain.set_boundary({'exterior': Br})
783
784        # Setup only one forcing term, time dependent rainfall
785        # restricted to a polygon enclosing triangle #1 (bce)
786        domain.forcing_terms = []
787        R = Rainfall(domain,
788                     rate=lambda t: 3*t + 7,
789                     polygon = [[1,1], [2,1], [2,2], [1,2]])
790
791        assert num.allclose(R.exchange_area, 2)
792       
793        domain.forcing_terms.append(R)
794
795        domain.time = 10.
796
797        domain.compute_forcing_terms()
798
799        assert num.allclose(domain.quantities['stage'].explicit_update[1],
800                            (3*domain.time + 7)/1000)
801        assert num.allclose(domain.quantities['stage'].explicit_update[0], 0)
802        assert num.allclose(domain.quantities['stage'].explicit_update[2:], 0)
803
804    def test_time_dependent_rainfall_using_starttime(self):
805        rainfall_poly = ensure_numeric([[1,1], [2,1], [2,2], [1,2]], num.float)
806
807        a = [0.0, 0.0]
808        b = [0.0, 2.0]
809        c = [2.0, 0.0]
810        d = [0.0, 4.0]
811        e = [2.0, 2.0]
812        f = [4.0, 0.0]
813
814        points = [a, b, c, d, e, f]
815        #             bac,     bce,     ecf,     dbe
816        vertices = [[1,0,2], [1,2,4], [4,2,5], [3,1,4]]
817
818        domain = Domain(points, vertices)
819
820        # Flat surface with 1m of water
821        domain.set_quantity('elevation', 0)
822        domain.set_quantity('stage', 1.0)
823        domain.set_quantity('friction', 0)
824
825        Br = Reflective_boundary(domain)
826        domain.set_boundary({'exterior': Br})
827
828        # Setup only one forcing term, time dependent rainfall
829        # restricted to a polygon enclosing triangle #1 (bce)
830        domain.forcing_terms = []
831        R = Rainfall(domain,
832                     rate=lambda t: 3*t + 7,
833                     polygon=rainfall_poly)                     
834
835        assert num.allclose(R.exchange_area, 2)
836       
837        domain.forcing_terms.append(R)
838
839        # This will test that time is set to starttime in set_starttime
840        domain.set_starttime(5.0)
841
842        domain.compute_forcing_terms()
843
844        assert num.allclose(domain.quantities['stage'].explicit_update[1],
845                            (3*domain.get_time() + 7)/1000)
846        assert num.allclose(domain.quantities['stage'].explicit_update[1],
847                            (3*domain.get_starttime() + 7)/1000)
848
849        assert num.allclose(domain.quantities['stage'].explicit_update[0], 0)
850        assert num.allclose(domain.quantities['stage'].explicit_update[2:], 0)
851
852    def test_time_dependent_rainfall_using_georef(self):
853        """test_time_dependent_rainfall_using_georef
854
855        This will also test the General forcing term using georef
856        """
857
858        # Mesh in zone 56 (absolute coords)
859        x0 = 314036.58727982
860        y0 = 6224951.2960092
861
862        rainfall_poly = ensure_numeric([[1,1], [2,1], [2,2], [1,2]], num.float)
863        rainfall_poly += [x0, y0]
864
865        a = [0.0, 0.0]
866        b = [0.0, 2.0]
867        c = [2.0, 0.0]
868        d = [0.0, 4.0]
869        e = [2.0, 2.0]
870        f = [4.0, 0.0]
871
872        points = [a, b, c, d, e, f]
873        #             bac,     bce,     ecf,     dbe
874        vertices = [[1,0,2], [1,2,4], [4,2,5], [3,1,4]]
875
876        domain = Domain(points, vertices,
877                        geo_reference=Geo_reference(56, x0, y0))
878
879        # Flat surface with 1m of water
880        domain.set_quantity('elevation', 0)
881        domain.set_quantity('stage', 1.0)
882        domain.set_quantity('friction', 0)
883
884        Br = Reflective_boundary(domain)
885        domain.set_boundary({'exterior': Br})
886
887        # Setup only one forcing term, time dependent rainfall
888        # restricted to a polygon enclosing triangle #1 (bce)
889        domain.forcing_terms = []
890        R = Rainfall(domain,
891                     rate=lambda t: 3*t + 7,
892                     polygon=rainfall_poly)
893
894        assert num.allclose(R.exchange_area, 2)
895       
896        domain.forcing_terms.append(R)
897
898        # This will test that time is set to starttime in set_starttime
899        domain.set_starttime(5.0)
900
901        domain.compute_forcing_terms()
902
903        assert num.allclose(domain.quantities['stage'].explicit_update[1],
904                            (3*domain.get_time() + 7)/1000)
905        assert num.allclose(domain.quantities['stage'].explicit_update[1],
906                            (3*domain.get_starttime() + 7)/1000)
907
908        assert num.allclose(domain.quantities['stage'].explicit_update[0], 0)
909        assert num.allclose(domain.quantities['stage'].explicit_update[2:], 0)
910
911    def test_time_dependent_rainfall_restricted_by_polygon_with_default(self):
912        """
913        Test that default rainfall can be used when given rate runs out of data.
914        """
915
916        import warnings
917        warnings.simplefilter('ignore', UserWarning)
918
919
920        a = [0.0, 0.0]
921        b = [0.0, 2.0]
922        c = [2.0, 0.0]
923        d = [0.0, 4.0]
924        e = [2.0, 2.0]
925        f = [4.0, 0.0]
926
927        points = [a, b, c, d, e, f]
928        #             bac,     bce,     ecf,     dbe
929        vertices = [[1,0,2], [1,2,4], [4,2,5], [3,1,4]]
930
931        domain = Domain(points, vertices)
932
933        # Flat surface with 1m of water
934        domain.set_quantity('elevation', 0)
935        domain.set_quantity('stage', 1.0)
936        domain.set_quantity('friction', 0)
937
938        Br = Reflective_boundary(domain)
939        domain.set_boundary({'exterior': Br})
940
941        # Setup only one forcing term, time dependent rainfall
942        # that expires at t==20
943        from anuga.fit_interpolate.interpolate import Modeltime_too_late
944
945        def main_rate(t):
946            if t > 20:
947                msg = 'Model time exceeded.'
948                raise Modeltime_too_late, msg
949            else:
950                return 3*t + 7
951
952        domain.forcing_terms = []
953        R = Rainfall(domain,
954                     rate=main_rate,
955                     polygon = [[1,1], [2,1], [2,2], [1,2]],
956                     default_rate=5.0)
957
958        assert num.allclose(R.exchange_area, 2)
959       
960        domain.forcing_terms.append(R)
961
962        domain.time = 10.
963
964        domain.compute_forcing_terms()
965
966        assert num.allclose(domain.quantities['stage'].explicit_update[1],
967                            (3*domain.time+7)/1000)
968        assert num.allclose(domain.quantities['stage'].explicit_update[0], 0)
969        assert num.allclose(domain.quantities['stage'].explicit_update[2:], 0)
970
971        domain.time = 100.
972        domain.quantities['stage'].explicit_update[:] = 0.0     # Reset
973        domain.compute_forcing_terms()
974
975        assert num.allclose(domain.quantities['stage'].explicit_update[1],
976                            5.0/1000) # Default value
977        assert num.allclose(domain.quantities['stage'].explicit_update[0], 0)
978        assert num.allclose(domain.quantities['stage'].explicit_update[2:], 0)
979
980    def test_rainfall_forcing_with_evolve(self):
981        """test_rainfall_forcing_with_evolve
982
983        Test how forcing terms are called within evolve
984        """
985
986        # FIXME(Ole): This test is just to experiment
987        import warnings
988        warnings.simplefilter('ignore', UserWarning)
989
990
991        a = [0.0, 0.0]
992        b = [0.0, 2.0]
993        c = [2.0, 0.0]
994        d = [0.0, 4.0]
995        e = [2.0, 2.0]
996        f = [4.0, 0.0]
997
998        points = [a, b, c, d, e, f]
999        #             bac,     bce,     ecf,     dbe
1000        vertices = [[1,0,2], [1,2,4], [4,2,5], [3,1,4]]
1001
1002        domain = Domain(points, vertices)
1003
1004        # Flat surface with 1m of water
1005        domain.set_quantity('elevation', 0)
1006        domain.set_quantity('stage', 1.0)
1007        domain.set_quantity('friction', 0)
1008
1009        Br = Reflective_boundary(domain)
1010        domain.set_boundary({'exterior': Br})
1011
1012        # Setup only one forcing term, time dependent rainfall
1013        # that expires at t==20
1014        from anuga.fit_interpolate.interpolate import Modeltime_too_late
1015
1016        def main_rate(t):
1017            if t > 20:
1018                msg = 'Model time exceeded.'
1019                raise Modeltime_too_late, msg
1020            else:
1021                return 3*t + 7
1022
1023        domain.forcing_terms = []
1024        R = Rainfall(domain,
1025                     rate=main_rate,
1026                     polygon=[[1,1], [2,1], [2,2], [1,2]],
1027                     default_rate=5.0)
1028
1029        assert num.allclose(R.exchange_area, 2)
1030       
1031        domain.forcing_terms.append(R)
1032
1033        for t in domain.evolve(yieldstep=1, finaltime=25):
1034            pass
1035            #FIXME(Ole):  A test here is hard because explicit_update also
1036            # receives updates from the flux calculation.
1037
1038
1039    def test_rainfall_forcing_with_evolve_1(self):
1040        """test_rainfall_forcing_with_evolve
1041
1042        Test how forcing terms are called within evolve.
1043        This test checks that proper exception is thrown when no default_rate is set
1044        """
1045
1046        import warnings
1047        warnings.simplefilter('ignore', UserWarning)
1048
1049
1050        a = [0.0, 0.0]
1051        b = [0.0, 2.0]
1052        c = [2.0, 0.0]
1053        d = [0.0, 4.0]
1054        e = [2.0, 2.0]
1055        f = [4.0, 0.0]
1056
1057        points = [a, b, c, d, e, f]
1058        #             bac,     bce,     ecf,     dbe
1059        vertices = [[1,0,2], [1,2,4], [4,2,5], [3,1,4]]
1060
1061        domain = Domain(points, vertices)
1062
1063        # Flat surface with 1m of water
1064        domain.set_quantity('elevation', 0)
1065        domain.set_quantity('stage', 1.0)
1066        domain.set_quantity('friction', 0)
1067
1068        Br = Reflective_boundary(domain)
1069        domain.set_boundary({'exterior': Br})
1070
1071        # Setup only one forcing term, time dependent rainfall
1072        # that expires at t==20
1073        from anuga.fit_interpolate.interpolate import Modeltime_too_late
1074
1075        def main_rate(t):
1076            if t > 20:
1077                msg = 'Model time exceeded.'
1078                raise Modeltime_too_late, msg
1079            else:
1080                return 3*t + 7
1081
1082        domain.forcing_terms = []
1083        R = Rainfall(domain,
1084                     rate=main_rate,
1085                     polygon=[[1,1], [2,1], [2,2], [1,2]])
1086
1087
1088        assert num.allclose(R.exchange_area, 2)
1089       
1090        domain.forcing_terms.append(R)
1091        #for t in domain.evolve(yieldstep=1, finaltime=25):
1092        #    pass
1093               
1094        try:
1095            for t in domain.evolve(yieldstep=1, finaltime=25):
1096                pass
1097        except Modeltime_too_late, e:
1098            # Test that error message is as expected
1099            assert 'can specify keyword argument default_rate in the forcing function' in str(e)
1100        else:
1101            raise Exception, 'Should have raised exception'
1102
1103    def test_constant_wind_stress_from_file(self):
1104        from anuga.config import rho_a, rho_w, eta_w
1105        from math import pi, cos, sin
1106
1107        cellsize = 25
1108        nrows=5; ncols = 6;
1109        refzone=50
1110        xllcorner=366000;yllcorner=6369500;
1111        number_of_timesteps = 6
1112        timestep=12*60
1113        eps=2e-16
1114
1115        points, vertices, boundary =rectangular(nrows-2,ncols-2,
1116                                                len1=cellsize*(ncols-1),
1117                                                len2=cellsize*(nrows-1),
1118                                                origin=(xllcorner,yllcorner))
1119
1120        domain = Domain(points, vertices, boundary)
1121        midpoints = domain.get_centroid_coordinates()
1122
1123        # Flat surface with 1m of water
1124        domain.set_quantity('elevation', 0)
1125        domain.set_quantity('stage', 1.0)
1126        domain.set_quantity('friction', 0)
1127
1128        Br = Reflective_boundary(domain)
1129        domain.set_boundary({'top': Br, 'bottom' :Br, 'left': Br, 'right': Br})
1130
1131        # Setup only one forcing term, constant wind stress
1132        s = 100
1133        phi = 135
1134        pressure=1000
1135        domain.forcing_terms = []
1136        field_sts_filename = 'wind_field'
1137        self.write_wind_pressure_field_sts(field_sts_filename,
1138                                      nrows=nrows,
1139                                      ncols=ncols,
1140                                      cellsize=cellsize,
1141                                      origin=(xllcorner,yllcorner),
1142                                      refzone=50,
1143                                      timestep=timestep,
1144                                      number_of_timesteps=10,
1145                                      speed=s,
1146                                      angle=phi,
1147                                      pressure=pressure)
1148
1149        sts2sww_mesh(field_sts_filename,spatial_thinning=1,
1150                     verbose=False)
1151
1152        # Setup wind stress
1153        F = file_function(field_sts_filename+'.sww', domain,
1154                          quantities=['wind_speed', 'wind_angle'],
1155                          interpolation_points = midpoints)
1156
1157        W = Wind_stress(F,use_coordinates=False)
1158        domain.forcing_terms.append(W)
1159        domain.compute_forcing_terms()
1160
1161        const = eta_w*rho_a / rho_w
1162
1163        # Convert to radians
1164        phi = phi*pi / 180
1165
1166        # Compute velocity vector (u, v)
1167        u = s*cos(phi)
1168        v = s*sin(phi)
1169
1170        # Compute wind stress
1171        S = const * num.sqrt(u**2 + v**2)
1172
1173        assert num.allclose(domain.quantities['stage'].explicit_update, 0)
1174        assert num.allclose(domain.quantities['xmomentum'].explicit_update, S*u)
1175        assert num.allclose(domain.quantities['ymomentum'].explicit_update, S*v)
1176
1177    def test_variable_windfield_from_file(self):
1178        from anuga.config import rho_a, rho_w, eta_w
1179        from math import pi, cos, sin
1180        from anuga.config import time_format
1181
1182        cellsize = 25
1183        #nrows=25; ncols = 25;
1184        nrows=10; ncols = 10;
1185        refzone=50
1186        xllcorner=366000;yllcorner=6369500;
1187        number_of_timesteps = 10
1188        timestep=1
1189        eps=2.e-16
1190        spatial_thinning=1
1191
1192        points, vertices, boundary =rectangular(nrows-2,ncols-2,
1193                                                len1=cellsize*(ncols-1),
1194                                                len2=cellsize*(nrows-1),
1195                                                origin=(xllcorner,yllcorner))
1196
1197        time=num.arange(0,10,1,num.float)
1198        eval_time=time[7];
1199
1200        domain = Domain(points, vertices, boundary)
1201        midpoints = domain.get_centroid_coordinates()
1202        vertexpoints = domain.get_nodes()
1203
1204        """
1205        x=grid_1d(xllcorner,cellsize,ncols)
1206        y=grid_1d(yllcorner,cellsize,nrows)
1207        X,Y=num.meshgrid(x,y)
1208        interpolation_points=num.empty((X.shape[0]*X.shape[1],2),num.float)
1209        k=0
1210        for i in range(X.shape[0]):
1211            for j in range(X.shape[1]):
1212                interpolation_points[k,0]=X[i,j]
1213                interpolation_points[k,1]=Y[i,j]
1214                k+=1
1215
1216        z=spatial_linear_varying_speed(eval_time,interpolation_points[:,0],
1217                                       interpolation_points[:,1])
1218
1219        k=0
1220        Z=num.empty((X.shape[0],X.shape[1]),num.float)
1221        for i in range(X.shape[0]):
1222            for j in range(X.shape[1]):
1223                Z[i,j]=z[k]
1224                k+=1
1225
1226        Q=num.empty((time.shape[0],points.shape[0]),num.float)
1227        for i, t in enumerate(time):
1228            Q[i,:]=spatial_linear_varying_speed(t,points[:,0],points[:,1])
1229
1230        from interpolate import Interpolation_function
1231        I  = Interpolation_function(time,Q,
1232                                    vertex_coordinates = points,
1233                                    triangles = domain.triangles,
1234                                    #interpolation_points = midpoints,
1235                                    interpolation_points=interpolation_points,
1236                                    verbose=False)
1237
1238        V=num.empty((X.shape[0],X.shape[1]),num.float)
1239        for k in range(len(interpolation_points)):
1240            assert num.allclose(I(eval_time,k),z[k])
1241            V[k/X.shape[1],k%X.shape[1]]=I(eval_time,k)
1242
1243
1244           import mpl_toolkits.mplot3d.axes3d as p3
1245           fig=P.figure()
1246           ax = p3.Axes3D(fig)
1247           ax.plot_surface(X,Y,V)
1248           ax.plot_surface(X,Y,Z)
1249           P.show()
1250
1251
1252        """
1253
1254        # Flat surface with 1m of water
1255        domain.set_quantity('elevation', 0)
1256        domain.set_quantity('stage', 1.0)
1257        domain.set_quantity('friction', 0)
1258
1259        domain.time = 7*timestep    # Take a time that is represented in file (not zero)
1260
1261        # Write wind stress file (ensure that domain.time is covered)
1262
1263        field_sts_filename = 'wind_field'
1264        self.write_wind_pressure_field_sts(field_sts_filename,
1265                                      nrows=nrows,
1266                                      ncols=ncols,
1267                                      cellsize=cellsize,
1268                                      origin=(xllcorner,yllcorner),
1269                                      refzone=50,
1270                                      timestep=timestep,
1271                                      number_of_timesteps=10,
1272                                      speed=spatial_linear_varying_speed,
1273                                      angle=spatial_linear_varying_angle,
1274                                      pressure=spatial_linear_varying_pressure)
1275
1276
1277        sts2sww_mesh(field_sts_filename,spatial_thinning=spatial_thinning,
1278                     verbose=False)
1279
1280        # Setup wind stress
1281        FW = file_function(field_sts_filename+'.sww', domain,
1282                          quantities=['wind_speed', 'wind_angle'],
1283                          interpolation_points = midpoints)
1284
1285        W = Wind_stress(FW,use_coordinates=False)
1286
1287        domain.forcing_terms = []
1288        domain.forcing_terms.append(W)
1289
1290        domain.compute_forcing_terms()
1291
1292        # Compute reference solution
1293        const = eta_w*rho_a / rho_w
1294
1295        N = len(domain)    # number_of_triangles
1296
1297        xc = domain.get_centroid_coordinates()
1298        t = domain.time
1299
1300        x = xc[:,0]
1301        y = xc[:,1]
1302        s_vec = spatial_linear_varying_speed(t,x,y)
1303        phi_vec = spatial_linear_varying_angle(t,x,y)
1304
1305        for k in range(N):
1306            # Convert to radians
1307            phi = phi_vec[k]*pi / 180
1308            s = s_vec[k]
1309
1310            # Compute velocity vector (u, v)
1311            u = s*cos(phi)
1312            v = s*sin(phi)
1313
1314            # Compute wind stress
1315            S = const * num.sqrt(u**2 + v**2)
1316
1317            assert num.allclose(domain.quantities['stage'].explicit_update[k],0)
1318
1319            assert num.allclose(domain.quantities['xmomentum'].\
1320                                    explicit_update[k],S*u,eps)
1321            assert num.allclose(domain.quantities['ymomentum'].\
1322                                     explicit_update[k],S*v,eps)
1323
1324        os.remove(field_sts_filename+'.sts')
1325        os.remove(field_sts_filename+'.sww')
1326
1327    def test_variable_pressurefield_from_file(self):
1328        from anuga.config import rho_a, rho_w, eta_w
1329        from math import pi, cos, sin
1330        from anuga.config import time_format
1331
1332        cellsize = 25
1333        #nrows=25; ncols = 25;
1334        nrows=10; ncols = 10;
1335        refzone=50
1336        xllcorner=366000;yllcorner=6369500;
1337        number_of_timesteps = 10
1338        timestep=1
1339        eps=2.e-16
1340        spatial_thinning=1
1341
1342        points, vertices, boundary =rectangular(nrows-2,ncols-2,
1343                                                len1=cellsize*(ncols-1),
1344                                                len2=cellsize*(nrows-1),
1345                                                origin=(xllcorner,yllcorner))
1346
1347        time=num.arange(0,10,1,num.float)
1348        eval_time=time[7];
1349
1350        domain = Domain(points, vertices, boundary)
1351        midpoints = domain.get_centroid_coordinates()
1352        vertexpoints = domain.get_nodes()
1353
1354        # Flat surface with 1m of water
1355        domain.set_quantity('elevation', 0)
1356        domain.set_quantity('stage', 1.0)
1357        domain.set_quantity('friction', 0)
1358
1359        domain.time = 7*timestep    # Take a time that is represented in file (not zero)
1360
1361        # Write wind stress file (ensure that domain.time is covered)
1362
1363        field_sts_filename = 'wind_field'
1364        self.write_wind_pressure_field_sts(field_sts_filename,
1365                                      nrows=nrows,
1366                                      ncols=ncols,
1367                                      cellsize=cellsize,
1368                                      origin=(xllcorner,yllcorner),
1369                                      refzone=50,
1370                                      timestep=timestep,
1371                                      number_of_timesteps=10,
1372                                      speed=spatial_linear_varying_speed,
1373                                      angle=spatial_linear_varying_angle,
1374                                      pressure=spatial_linear_varying_pressure)
1375
1376
1377        sts2sww_mesh(field_sts_filename,spatial_thinning=spatial_thinning,
1378                     verbose=False)
1379
1380        # Setup barometric pressure
1381        FP = file_function(field_sts_filename+'.sww', domain,
1382                           quantities=['barometric_pressure'],
1383                           interpolation_points = vertexpoints)
1384
1385        P = Barometric_pressure(FP,use_coordinates=False)
1386
1387
1388        domain.forcing_terms = []
1389        domain.forcing_terms.append(P)
1390
1391        domain.compute_forcing_terms()
1392
1393        N = len(domain)    # number_of_triangles
1394
1395        xc = domain.get_centroid_coordinates()
1396        t = domain.time
1397
1398        x = xc[:,0]
1399        y = xc[:,1]
1400        p_vec = spatial_linear_varying_pressure(t,x,y)
1401
1402        h=1 #depth
1403        px=0.000025  #pressure gradient in x-direction
1404        py=0.0000125 #pressure gradient in y-direction
1405        for k in range(N):
1406            # Convert to radians
1407            p = p_vec[k]
1408
1409            assert num.allclose(domain.quantities['stage'].explicit_update[k],0)
1410
1411            assert num.allclose(domain.quantities['xmomentum'].\
1412                                    explicit_update[k],h*px/rho_w)
1413
1414            assert num.allclose(domain.quantities['ymomentum'].\
1415                                     explicit_update[k],h*py/rho_w)
1416
1417        os.remove(field_sts_filename+'.sts')
1418        os.remove(field_sts_filename+'.sww')
1419
1420    def test_constant_wind_stress_from_file_evolve(self):
1421        from anuga.config import rho_a, rho_w, eta_w
1422        from math import pi, cos, sin
1423        from anuga.config import time_format
1424
1425        cellsize = 25
1426        nrows=5; ncols = 6;
1427        refzone=50
1428        xllcorner=366000;yllcorner=6369500;
1429        number_of_timesteps = 27
1430        timestep=1
1431        eps=2e-16
1432
1433        points, vertices, boundary =rectangular(nrows-2,ncols-2,
1434                                                len1=cellsize*(ncols-1),
1435                                                len2=cellsize*(nrows-1),
1436                                                origin=(xllcorner,yllcorner))
1437
1438        domain = Domain(points, vertices, boundary)
1439        midpoints = domain.get_centroid_coordinates()
1440
1441        # Flat surface with 1m of water
1442        domain.set_quantity('elevation', 0)
1443        domain.set_quantity('stage', 1.0)
1444        domain.set_quantity('friction', 0)
1445
1446        Br = Reflective_boundary(domain)
1447        domain.set_boundary({'top': Br, 'bottom' :Br, 'left': Br, 'right': Br})
1448
1449        # Setup only one forcing term, constant wind stress
1450        s = 100
1451        phi = 135
1452        field_sts_filename = 'wind_field'
1453        self.write_wind_pressure_field_sts(field_sts_filename,
1454                                      nrows=nrows,
1455                                      ncols=ncols,
1456                                      cellsize=cellsize,
1457                                      origin=(xllcorner,yllcorner),
1458                                      refzone=50,
1459                                      timestep=timestep,
1460                                      number_of_timesteps=number_of_timesteps,
1461                                      speed=s,
1462                                      angle=phi)
1463
1464        sts2sww_mesh(field_sts_filename,spatial_thinning=1,
1465                     verbose=False)
1466
1467        # Setup wind stress
1468        F = file_function(field_sts_filename+'.sww', domain,
1469                          quantities=['wind_speed', 'wind_angle'],
1470                          interpolation_points = midpoints)
1471
1472        W = Wind_stress(F,use_coordinates=False)
1473        domain.forcing_terms.append(W)
1474
1475        valuesUsingFunction=num.empty((3,number_of_timesteps+1,midpoints.shape[0]),
1476                                      num.float)
1477        i=0
1478        for t in domain.evolve(yieldstep=1, finaltime=number_of_timesteps*timestep):
1479            valuesUsingFunction[0,i]=domain.quantities['stage'].explicit_update
1480            valuesUsingFunction[1,i]=domain.quantities['xmomentum'].explicit_update
1481            valuesUsingFunction[2,i]=domain.quantities['ymomentum'].explicit_update
1482            i+=1
1483
1484
1485        domain_II = Domain(points, vertices, boundary)
1486
1487        # Flat surface with 1m of water
1488        domain_II.set_quantity('elevation', 0)
1489        domain_II.set_quantity('stage', 1.0)
1490        domain_II.set_quantity('friction', 0)
1491
1492        Br = Reflective_boundary(domain_II)
1493        domain_II.set_boundary({'top': Br, 'bottom' :Br, 'left': Br, 'right': Br})
1494
1495        s = 100
1496        phi = 135
1497        domain_II.forcing_terms = []
1498        domain_II.forcing_terms.append(Wind_stress(s, phi))
1499
1500        i=0;
1501        for t in domain_II.evolve(yieldstep=1, 
1502                                  finaltime=number_of_timesteps*timestep):
1503            assert num.allclose(valuesUsingFunction[0,i],domain_II.quantities['stage'].explicit_update), max(valuesUsingFunction[0,i]-domain_II.quantities['stage'].explicit_update)
1504            assert  num.allclose(valuesUsingFunction[1,i],domain_II.quantities['xmomentum'].explicit_update)
1505            assert num.allclose(valuesUsingFunction[2,i],domain_II.quantities['ymomentum'].explicit_update)
1506            i+=1
1507
1508        os.remove(field_sts_filename+'.sts')
1509        os.remove(field_sts_filename+'.sww')
1510
1511    def test_temporally_varying_wind_stress_from_file_evolve(self):
1512        from anuga.config import rho_a, rho_w, eta_w
1513        from math import pi, cos, sin
1514        from anuga.config import time_format
1515
1516        cellsize = 25
1517        #nrows=20; ncols = 20;
1518        nrows=10; ncols = 10;
1519        refzone=50
1520        xllcorner=366000;yllcorner=6369500;
1521        number_of_timesteps = 28
1522        timestep=1.
1523        eps=2e-16
1524
1525        #points, vertices, boundary =rectangular(10,10,
1526        points, vertices, boundary =rectangular(5,5,
1527                                                len1=cellsize*(ncols-1),
1528                                                len2=cellsize*(nrows-1),
1529                                                origin=(xllcorner,yllcorner))
1530
1531        domain = Domain(points, vertices, boundary)
1532        midpoints = domain.get_centroid_coordinates()
1533
1534        # Flat surface with 1m of water
1535        domain.set_quantity('elevation', 0)
1536        domain.set_quantity('stage', 1.0)
1537        domain.set_quantity('friction', 0)
1538
1539        Br = Reflective_boundary(domain)
1540        domain.set_boundary({'top': Br, 'bottom' :Br, 'left': Br, 'right': Br})
1541
1542        # Setup only one forcing term, constant wind stress
1543        field_sts_filename = 'wind_field'
1544        self.write_wind_pressure_field_sts(field_sts_filename,
1545                                      nrows=nrows,
1546                                      ncols=ncols,
1547                                      cellsize=cellsize,
1548                                      origin=(xllcorner,yllcorner),
1549                                      refzone=50,
1550                                      timestep=timestep,
1551                                      number_of_timesteps=number_of_timesteps,
1552                                      speed=time_varying_speed,
1553                                      angle=time_varying_angle,
1554                                      pressure=time_varying_pressure)
1555
1556        sts2sww_mesh(field_sts_filename,spatial_thinning=1,
1557                     verbose=False)
1558
1559        # Setup wind stress
1560        F = file_function(field_sts_filename+'.sww', domain,
1561                          quantities=['wind_speed', 'wind_angle'],
1562                          interpolation_points = midpoints)
1563
1564        #W = Wind_stress(F,use_coordinates=False)
1565        W = Wind_stress_fast(F,filename=field_sts_filename+'.sww', domain=domain)
1566        domain.forcing_terms.append(W)
1567
1568        valuesUsingFunction=num.empty((3,2*number_of_timesteps,midpoints.shape[0]),
1569                                      num.float)
1570        i=0
1571        for t in domain.evolve(yieldstep=timestep/2., finaltime=(number_of_timesteps-1)*timestep):
1572            valuesUsingFunction[0,i]=domain.quantities['stage'].explicit_update
1573            valuesUsingFunction[1,i]=domain.quantities['xmomentum'].explicit_update
1574            valuesUsingFunction[2,i]=domain.quantities['ymomentum'].explicit_update
1575            i+=1
1576
1577
1578        domain_II = Domain(points, vertices, boundary)
1579
1580        # Flat surface with 1m of water
1581        domain_II.set_quantity('elevation', 0)
1582        domain_II.set_quantity('stage', 1.0)
1583        domain_II.set_quantity('friction', 0)
1584
1585        Br = Reflective_boundary(domain_II)
1586        domain_II.set_boundary({'top': Br, 'bottom' :Br, 'left': Br, 'right': Br})
1587
1588        domain_II.forcing_terms.append(Wind_stress(s=time_varying_speed, 
1589                                                   phi=time_varying_angle))
1590
1591        i=0;
1592        for t in domain_II.evolve(yieldstep=timestep/2., 
1593                                  finaltime=(number_of_timesteps-1)*timestep):
1594            assert num.allclose(valuesUsingFunction[0,i],
1595                                domain_II.quantities['stage'].explicit_update,
1596                                eps)
1597            #print i,valuesUsingFunction[1,i]
1598            assert  num.allclose(valuesUsingFunction[1,i],
1599                                 domain_II.quantities['xmomentum'].explicit_update,
1600                                 eps),(valuesUsingFunction[1,i]-
1601                                 domain_II.quantities['xmomentum'].explicit_update)
1602            assert num.allclose(valuesUsingFunction[2,i],
1603                                domain_II.quantities['ymomentum'].explicit_update,
1604                                eps)
1605            #if i==1: assert-1==1
1606            i+=1
1607
1608        os.remove(field_sts_filename+'.sts')
1609        os.remove(field_sts_filename+'.sww')
1610
1611    def test_spatially_varying_wind_stress_from_file_evolve(self):
1612        from anuga.config import rho_a, rho_w, eta_w
1613        from math import pi, cos, sin
1614        from anuga.config import time_format
1615
1616        cellsize = 25
1617        nrows=20; ncols = 20;
1618        nrows=10; ncols = 10;
1619        refzone=50
1620        xllcorner=366000;yllcorner=6369500;
1621        number_of_timesteps = 28
1622        timestep=1.
1623        eps=2e-16
1624
1625        #points, vertices, boundary =rectangular(10,10,
1626        points, vertices, boundary =rectangular(5,5,
1627                                                len1=cellsize*(ncols-1),
1628                                                len2=cellsize*(nrows-1),
1629                                                origin=(xllcorner,yllcorner))
1630
1631        domain = Domain(points, vertices, boundary)
1632        midpoints = domain.get_centroid_coordinates()
1633
1634        # Flat surface with 1m of water
1635        domain.set_quantity('elevation', 0)
1636        domain.set_quantity('stage', 1.0)
1637        domain.set_quantity('friction', 0)
1638
1639        Br = Reflective_boundary(domain)
1640        domain.set_boundary({'top': Br, 'bottom' :Br, 'left': Br, 'right': Br})
1641
1642        # Setup only one forcing term, constant wind stress
1643        field_sts_filename = 'wind_field'
1644        self.write_wind_pressure_field_sts(field_sts_filename,
1645                                      nrows=nrows,
1646                                      ncols=ncols,
1647                                      cellsize=cellsize,
1648                                      origin=(xllcorner,yllcorner),
1649                                      refzone=50,
1650                                      timestep=timestep,
1651                                      number_of_timesteps=number_of_timesteps,
1652                                      speed=spatial_linear_varying_speed,
1653                                      angle=spatial_linear_varying_angle,
1654                                      pressure=spatial_linear_varying_pressure)
1655
1656        sts2sww_mesh(field_sts_filename,spatial_thinning=1,
1657                     verbose=False)
1658
1659        # Setup wind stress
1660        F = file_function(field_sts_filename+'.sww', domain,
1661                          quantities=['wind_speed', 'wind_angle'],
1662                          interpolation_points = midpoints)
1663
1664        W = Wind_stress(F,use_coordinates=False)
1665        domain.forcing_terms.append(W)
1666
1667        valuesUsingFunction=num.empty((3,number_of_timesteps,midpoints.shape[0]),
1668                                      num.float)
1669        i=0
1670        for t in domain.evolve(yieldstep=timestep, finaltime=(number_of_timesteps-1)*timestep):
1671            valuesUsingFunction[0,i]=domain.quantities['stage'].explicit_update
1672            valuesUsingFunction[1,i]=domain.quantities['xmomentum'].explicit_update
1673            valuesUsingFunction[2,i]=domain.quantities['ymomentum'].explicit_update
1674            i+=1
1675
1676
1677        domain_II = Domain(points, vertices, boundary)
1678
1679        # Flat surface with 1m of water
1680        domain_II.set_quantity('elevation', 0)
1681        domain_II.set_quantity('stage', 1.0)
1682        domain_II.set_quantity('friction', 0)
1683
1684        Br = Reflective_boundary(domain_II)
1685        domain_II.set_boundary({'top': Br, 'bottom' :Br, 'left': Br, 'right': Br})
1686
1687        domain_II.forcing_terms.append(Wind_stress(s=spatial_linear_varying_speed, 
1688                                                   phi=spatial_linear_varying_angle))
1689
1690        i=0;
1691        for t in domain_II.evolve(yieldstep=timestep, 
1692                                  finaltime=(number_of_timesteps-1)*timestep):
1693            #print valuesUsingFunction[1,i],domain_II.quantities['xmomentum'].explicit_update
1694            assert num.allclose(valuesUsingFunction[0,i],
1695                                domain_II.quantities['stage'].explicit_update,
1696                                eps)
1697            assert  num.allclose(valuesUsingFunction[1,i],
1698                                 domain_II.quantities['xmomentum'].explicit_update,
1699                                 eps)
1700            assert num.allclose(valuesUsingFunction[2,i],
1701                                domain_II.quantities['ymomentum'].explicit_update,
1702                                eps)
1703            i+=1
1704
1705        os.remove(field_sts_filename+'.sts')
1706        os.remove(field_sts_filename+'.sww')
1707
1708    def test_temporally_varying_pressure_stress_from_file_evolve(self):
1709        from anuga.config import rho_a, rho_w, eta_w
1710        from math import pi, cos, sin
1711        from anuga.config import time_format
1712
1713        cellsize = 25
1714        #nrows=20; ncols = 20;
1715        nrows=10; ncols = 10;
1716        refzone=50
1717        xllcorner=366000;yllcorner=6369500;
1718        number_of_timesteps = 28
1719        timestep=10.
1720        eps=2e-16
1721
1722        #print "Building mesh"
1723        #points, vertices, boundary =rectangular(10,10,
1724        points, vertices, boundary =rectangular(5,5,
1725                                                len1=cellsize*(ncols-1),
1726                                                len2=cellsize*(nrows-1),
1727                                                origin=(xllcorner,yllcorner))
1728
1729        domain = Domain(points, vertices, boundary)
1730        vertexpoints = domain.get_nodes()
1731
1732        # Flat surface with 1m of water
1733        domain.set_quantity('elevation', 0)
1734        domain.set_quantity('stage', 1.0)
1735        domain.set_quantity('friction', 0)
1736
1737        Br = Reflective_boundary(domain)
1738        domain.set_boundary({'top': Br, 'bottom' :Br, 'left': Br, 'right': Br})
1739
1740        # Setup only one forcing term, constant wind stress
1741        field_sts_filename = 'wind_field'
1742        #print 'Writing pressure field sts file'
1743        self.write_wind_pressure_field_sts(field_sts_filename,
1744                                      nrows=nrows,
1745                                      ncols=ncols,
1746                                      cellsize=cellsize,
1747                                      origin=(xllcorner,yllcorner),
1748                                      refzone=50,
1749                                      timestep=timestep,
1750                                      number_of_timesteps=number_of_timesteps,
1751                                      speed=time_varying_speed,
1752                                      angle=time_varying_angle,
1753                                      pressure=time_varying_pressure)
1754
1755        #print "converting sts to sww"
1756        sts2sww_mesh(field_sts_filename,spatial_thinning=1,
1757                     verbose=False)
1758
1759        #print 'initialising file_function'
1760        # Setup wind stress
1761        F = file_function(field_sts_filename+'.sww', domain,
1762                          quantities=['barometric_pressure'],
1763                          interpolation_points = vertexpoints)
1764
1765        #P = Barometric_pressure(F,use_coordinates=False)
1766        #print 'initialising pressure forcing term'
1767        P = Barometric_pressure_fast(p=F,filename=field_sts_filename+'.sww',domain=domain)
1768        domain.forcing_terms.append(P)
1769
1770        valuesUsingFunction=num.empty((3,2*number_of_timesteps,len(domain)),
1771                                      num.float)
1772        i=0
1773        import time as timer
1774        t0=timer.time()
1775        for t in domain.evolve(yieldstep=timestep/2., finaltime=(number_of_timesteps-1)*timestep):
1776            valuesUsingFunction[0,i]=domain.quantities['stage'].explicit_update
1777            valuesUsingFunction[1,i]=domain.quantities['xmomentum'].explicit_update
1778            valuesUsingFunction[2,i]=domain.quantities['ymomentum'].explicit_update
1779            i+=1
1780            #domain.write_time()
1781        t1=timer.time()
1782        #print "That took %fs seconds" %(t1-t0)
1783
1784
1785        domain_II = Domain(points, vertices, boundary)
1786
1787        # Flat surface with 1m of water
1788        domain_II.set_quantity('elevation', 0)
1789        domain_II.set_quantity('stage', 1.0)
1790        domain_II.set_quantity('friction', 0)
1791
1792        Br = Reflective_boundary(domain_II)
1793        domain_II.set_boundary({'top': Br, 'bottom' :Br, 'left': Br, 'right': Br})
1794
1795        domain_II.forcing_terms.append(Barometric_pressure(p=time_varying_pressure))
1796
1797        i=0;
1798        for t in domain_II.evolve(yieldstep=timestep/2., 
1799                                  finaltime=(number_of_timesteps-1)*timestep):
1800            assert num.allclose(valuesUsingFunction[0,i],
1801                                domain_II.quantities['stage'].explicit_update,
1802                                eps)
1803            assert  num.allclose(valuesUsingFunction[1,i],
1804                                 domain_II.quantities['xmomentum'].explicit_update,
1805                                 eps)
1806            assert num.allclose(valuesUsingFunction[2,i],
1807                                domain_II.quantities['ymomentum'].explicit_update,
1808                                eps)
1809            i+=1
1810
1811        os.remove(field_sts_filename+'.sts')
1812        os.remove(field_sts_filename+'.sww')
1813
1814    def test_spatially_varying_pressure_stress_from_file_evolve(self):
1815        from anuga.config import rho_a, rho_w, eta_w
1816        from math import pi, cos, sin
1817        from anuga.config import time_format
1818
1819        cellsize = 25
1820        #nrows=20; ncols = 20;
1821        nrows=10; ncols = 10;
1822        refzone=50
1823        xllcorner=366000;yllcorner=6369500;
1824        number_of_timesteps = 28
1825        timestep=1.
1826        eps=2e-16
1827
1828        #points, vertices, boundary =rectangular(10,10,
1829        points, vertices, boundary =rectangular(5,5,
1830                                                len1=cellsize*(ncols-1),
1831                                                len2=cellsize*(nrows-1),
1832                                                origin=(xllcorner,yllcorner))
1833
1834        domain = Domain(points, vertices, boundary)
1835        vertexpoints = domain.get_nodes()
1836
1837        # Flat surface with 1m of water
1838        domain.set_quantity('elevation', 0)
1839        domain.set_quantity('stage', 1.0)
1840        domain.set_quantity('friction', 0)
1841
1842        Br = Reflective_boundary(domain)
1843        domain.set_boundary({'top': Br, 'bottom' :Br, 'left': Br, 'right': Br})
1844
1845        # Setup only one forcing term, constant wind stress
1846        field_sts_filename = 'wind_field'
1847        self.write_wind_pressure_field_sts(field_sts_filename,
1848                                      nrows=nrows,
1849                                      ncols=ncols,
1850                                      cellsize=cellsize,
1851                                      origin=(xllcorner,yllcorner),
1852                                      refzone=50,
1853                                      timestep=timestep,
1854                                      number_of_timesteps=number_of_timesteps,
1855                                      speed=spatial_linear_varying_speed,
1856                                      angle=spatial_linear_varying_angle,
1857                                      pressure=spatial_linear_varying_pressure)
1858
1859        sts2sww_mesh(field_sts_filename,spatial_thinning=1,
1860                     verbose=False)
1861
1862        # Setup wind stress
1863        F = file_function(field_sts_filename+'.sww', domain,
1864                          quantities=['barometric_pressure'],
1865                          interpolation_points = vertexpoints)
1866
1867        P = Barometric_pressure(F,use_coordinates=False)
1868        domain.forcing_terms.append(P)
1869
1870        valuesUsingFunction=num.empty((3,number_of_timesteps,len(domain)),
1871                                      num.float)
1872        i=0
1873        for t in domain.evolve(yieldstep=timestep, finaltime=(number_of_timesteps-1)*timestep):
1874            valuesUsingFunction[0,i]=domain.quantities['stage'].explicit_update
1875            valuesUsingFunction[1,i]=domain.quantities['xmomentum'].explicit_update
1876            valuesUsingFunction[2,i]=domain.quantities['ymomentum'].explicit_update
1877            i+=1
1878
1879
1880        domain_II = Domain(points, vertices, boundary)
1881
1882        # Flat surface with 1m of water
1883        domain_II.set_quantity('elevation', 0)
1884        domain_II.set_quantity('stage', 1.0)
1885        domain_II.set_quantity('friction', 0)
1886
1887        Br = Reflective_boundary(domain_II)
1888        domain_II.set_boundary({'top': Br, 'bottom' :Br, 'left': Br, 'right': Br})
1889
1890        domain_II.forcing_terms.append(Barometric_pressure(p=spatial_linear_varying_pressure))
1891
1892        i=0;
1893        for t in domain_II.evolve(yieldstep=timestep, 
1894                                  finaltime=(number_of_timesteps-1)*timestep):
1895
1896            assert num.allclose(valuesUsingFunction[0,i],
1897                                domain_II.quantities['stage'].explicit_update,
1898                                eps)
1899            assert  num.allclose(valuesUsingFunction[1,i],
1900                                 domain_II.quantities['xmomentum'].explicit_update,
1901                                 eps)
1902            assert num.allclose(valuesUsingFunction[2,i],
1903                                domain_II.quantities['ymomentum'].explicit_update,
1904                                eps)
1905            i+=1
1906
1907        os.remove(field_sts_filename+'.sts')
1908        os.remove(field_sts_filename+'.sww')
1909
1910
1911           
1912if __name__ == "__main__":
1913    suite = unittest.makeSuite(Test_Forcing, 'test')
1914    runner = unittest.TextTestRunner(verbosity=1)
1915    runner.run(suite)
Note: See TracBrowser for help on using the repository browser.