# source:anuga_core/source/anuga/culvert_flows/test_culvert_class.py@7562

Last change on this file since 7562 was 7562, checked in by steve, 14 years ago

Updating the balanced and parallel code

File size: 30.0 KB
Line
1#!/usr/bin/env python
2
3
4import unittest
5import os.path
6import sys
7
8from anuga.utilities.system_tools import get_pathname_from_package
9from anuga.utilities.polygon import Polygon_function
10
11from anuga.abstract_2d_finite_volumes.mesh_factory import rectangular_cross
12from anuga.abstract_2d_finite_volumes.quantity import Quantity
13
14from anuga.shallow_water import Domain, Reflective_boundary,\
15    Dirichlet_boundary,\
16    Transmissive_boundary, Time_boundary
17
18from anuga.culvert_flows.culvert_class import Culvert_flow, Culvert_flow_rating, Culvert_flow_energy
19from anuga.culvert_flows.culvert_routines import boyd_generalised_culvert_model
20
21from math import pi,pow,sqrt
22
23import numpy as num
24
25
26# Helper functions
27def run_culvert_flow_problem(depth):
28    """Run flow with culvert given depth
29    """
30
31    length = 40.
32    width = 5.
33
34    dx = dy = 1           # Resolution: Length of subdivisions on both axes
35
36    points, vertices, boundary = rectangular_cross(int(length/dx),
37                                                   int(width/dy),
38                                                   len1=length,
39                                                   len2=width)
40    domain = Domain(points, vertices, boundary)
41    domain.set_name('Test_culvert_shallow') # Output name
42    domain.set_default_order(2)
43
44
45    #----------------------------------------------------------------------
46    # Setup initial conditions
47    #----------------------------------------------------------------------
48
49    def topography(x, y):
50        """Set up a weir
51
52        A culvert will connect either side
53        """
54        # General Slope of Topography
55        z=-x/1000
56
57        N = len(x)
58        for i in range(N):
59
60           # Sloping Embankment Across Channel
61            if 5.0 < x[i] < 10.1:
62                # Cut Out Segment for Culvert face
63                if  1.0+(x[i]-5.0)/5.0 <  y[i]  < 4.0 - (x[i]-5.0)/5.0:
64                   z[i]=z[i]
65                else:
66                   z[i] +=  0.5*(x[i] -5.0)    # Sloping Segment  U/S Face
67            if 10.0 < x[i] < 12.1:
68               z[i] +=  2.5                    # Flat Crest of Embankment
69            if 12.0 < x[i] < 14.5:
70                # Cut Out Segment for Culvert face
71                if  2.0-(x[i]-12.0)/2.5 <  y[i]  < 3.0 + (x[i]-12.0)/2.5:
72                   z[i]=z[i]
73                else:
74                   z[i] +=  2.5-1.0*(x[i] -12.0) # Sloping D/S Face
75
76
77        return z
78
79
80    domain.set_quantity('elevation', topography)
81    domain.set_quantity('friction', 0.01)         # Constant friction
82    domain.set_quantity('stage',
83                        expression='elevation + %f' % depth) # Shallow initial condition
84
85    # Boyd culvert
86    culvert = Culvert_flow(domain,
87                           label='Culvert No. 1',
88                           description='This culvert is a test unit 1.2m Wide by 0.75m High',
89                           end_point0=[9.0, 2.5],
90                           end_point1=[13.0, 2.5],
91                           width=1.20, height=0.75,
92                           culvert_routine=boyd_generalised_culvert_model,
93                           number_of_barrels=1,
94                           update_interval=2,
95                           verbose=False)
96
97
98    domain.forcing_terms.append(culvert)
99
100
101    #-----------------------------------------------------------------------
102    # Setup boundary conditions
103    #-----------------------------------------------------------------------
104
105    # Inflow based on Flow Depth and Approaching Momentum
106
107    Br = Reflective_boundary(domain)              # Solid reflective wall
108    domain.set_boundary({'left': Br, 'right': Br, 'top': Br, 'bottom': Br})
109
110
111
112    #-----------------------------------------------------------------------
113    # Evolve system through time
114    #-----------------------------------------------------------------------
115
116    #print 'depth', depth
117    ref_volume = domain.get_quantity('stage').get_integral()
118    for t in domain.evolve(yieldstep = 0.1, finaltime = 25):
119        new_volume = domain.get_quantity('stage').get_integral()
120
121        msg = ('Total volume has changed: Is %.8f m^3 should have been %.8f m^3'
122               % (new_volume, ref_volume))
123        assert num.allclose(new_volume, ref_volume), msg
124
125
126    os.remove('Test_culvert_shallow.sww')
127
128class Test_Culvert(unittest.TestCase):
129    def setUp(self):
130        pass
131
132    def tearDown(self):
133        pass
134
135
136    def test_that_culvert_runs_rating(self):
137        """test_that_culvert_runs_rating
138
139        This test exercises the culvert and checks values outside rating curve
140        are dealt with
141        """
142
143        path = get_pathname_from_package('anuga.culvert_flows')
144
145        length = 40.
146        width = 5.
147
148        dx = dy = 1           # Resolution: Length of subdivisions on both axes
149
150        points, vertices, boundary = rectangular_cross(int(length/dx),
151                                                       int(width/dy),
152                                                       len1=length,
153                                                       len2=width)
154        domain = Domain(points, vertices, boundary)
155        domain.set_name('Test_culvert')                 # Output name
156        domain.set_default_order(2)
157
158
159        #----------------------------------------------------------------------
160        # Setup initial conditions
161        #----------------------------------------------------------------------
162
163        def topography(x, y):
164            """Set up a weir
165
166            A culvert will connect either side
167            """
168            # General Slope of Topography
169            z=-x/1000
170
171            N = len(x)
172            for i in range(N):
173
174               # Sloping Embankment Across Channel
175                if 5.0 < x[i] < 10.1:
176                    # Cut Out Segment for Culvert face
177                    if  1.0+(x[i]-5.0)/5.0 <  y[i]  < 4.0 - (x[i]-5.0)/5.0:
178                       z[i]=z[i]
179                    else:
180                       z[i] +=  0.5*(x[i] -5.0)    # Sloping Segment  U/S Face
181                if 10.0 < x[i] < 12.1:
182                   z[i] +=  2.5                    # Flat Crest of Embankment
183                if 12.0 < x[i] < 14.5:
184                    # Cut Out Segment for Culvert face
185                    if  2.0-(x[i]-12.0)/2.5 <  y[i]  < 3.0 + (x[i]-12.0)/2.5:
186                       z[i]=z[i]
187                    else:
188                       z[i] +=  2.5-1.0*(x[i] -12.0) # Sloping D/S Face
189
190
191            return z
192
193
194        domain.set_quantity('elevation', topography)
195        domain.set_quantity('friction', 0.01)         # Constant friction
196        domain.set_quantity('stage',
197                            expression='elevation')   # Dry initial condition
198
199        filename=os.path.join(path, 'example_rating_curve.csv')
200        culvert = Culvert_flow(domain,
201                               culvert_description_filename=filename,
202                               end_point0=[9.0, 2.5],
203                               end_point1=[13.0, 2.5],
204                               width=1.00,
206                               verbose=False)
207
208        domain.forcing_terms.append(culvert)
209
210
211        #-----------------------------------------------------------------------
212        # Setup boundary conditions
213        #-----------------------------------------------------------------------
214
215        # Inflow based on Flow Depth and Approaching Momentum
216        Bi = Dirichlet_boundary([0.0, 0.0, 0.0])
217        Br = Reflective_boundary(domain)              # Solid reflective wall
218        Bo = Dirichlet_boundary([-5, 0, 0])           # Outflow
219
220        # Upstream and downstream conditions that will exceed the rating curve
221        # I.e produce delta_h outside the range [0, 10] specified in the the
222        # file example_rating_curve.csv
223        Btus = Time_boundary(domain, lambda t: [100*num.sin(2*pi*(t-4)/10), 0.0, 0.0])
224        Btds = Time_boundary(domain, lambda t: [-5*(num.cos(2*pi*(t-4)/20)), 0.0, 0.0])
225        domain.set_boundary({'left': Btus, 'right': Btds, 'top': Br, 'bottom': Br})
226
227
228        #-----------------------------------------------------------------------
229        # Evolve system through time
230        #-----------------------------------------------------------------------
231
232        min_delta_w = sys.maxint
233        max_delta_w = -min_delta_w
234        for t in domain.evolve(yieldstep = 1, finaltime = 25):
235            delta_w = culvert.inlet.stage - culvert.outlet.stage
236
237            if delta_w > max_delta_w: max_delta_w = delta_w
238            if delta_w < min_delta_w: min_delta_w = delta_w
239
240            pass
241
242        # Check that extreme values in rating curve have been exceeded
243        # so that we know that condition has been exercised
244        assert min_delta_w < 0
245        assert max_delta_w > 10
246
247
248        os.remove('Test_culvert.sww')
249
250    def test_that_culvert_dry_bed_rating_does_not_produce_flow(self):
251        """test_that_culvert_in_dry_bed_does_not_produce_flow(self):
252
253        Test that culvert on a sloping dry bed doesn't produce flows
254        although there will be a 'pressure' head due to delta_w > 0
255
256        This one is using the rating curve variant
257        """
258
259        path = get_pathname_from_package('anuga.culvert_flows')
260
261        length = 40.
262        width = 5.
263
264        dx = dy = 1           # Resolution: Length of subdivisions on both axes
265
266        points, vertices, boundary = rectangular_cross(int(length/dx),
267                                                       int(width/dy),
268                                                       len1=length,
269                                                       len2=width)
270        domain = Domain(points, vertices, boundary)
271        domain.set_name('Test_culvert_dry') # Output name
272        domain.set_default_order(2)
273
274
275        #----------------------------------------------------------------------
276        # Setup initial conditions
277        #----------------------------------------------------------------------
278
279        def topography(x, y):
280            """Set up a weir
281
282            A culvert will connect either side
283            """
284            # General Slope of Topography
285            z=-x/1000
286
287            N = len(x)
288            for i in range(N):
289
290               # Sloping Embankment Across Channel
291                if 5.0 < x[i] < 10.1:
292                    # Cut Out Segment for Culvert face
293                    if  1.0+(x[i]-5.0)/5.0 <  y[i]  < 4.0 - (x[i]-5.0)/5.0:
294                       z[i]=z[i]
295                    else:
296                       z[i] +=  0.5*(x[i] -5.0)    # Sloping Segment  U/S Face
297                if 10.0 < x[i] < 12.1:
298                   z[i] +=  2.5                    # Flat Crest of Embankment
299                if 12.0 < x[i] < 14.5:
300                    # Cut Out Segment for Culvert face
301                    if  2.0-(x[i]-12.0)/2.5 <  y[i]  < 3.0 + (x[i]-12.0)/2.5:
302                       z[i]=z[i]
303                    else:
304                       z[i] +=  2.5-1.0*(x[i] -12.0) # Sloping D/S Face
305
306
307            return z
308
309
310        domain.set_quantity('elevation', topography)
311        domain.set_quantity('friction', 0.01)         # Constant friction
312        domain.set_quantity('stage',
313                            expression='elevation')   # Dry initial condition
314
315
316        filename = os.path.join(path, 'example_rating_curve.csv')
317        culvert = Culvert_flow(domain,
318                               culvert_description_filename=filename,
319                               end_point0=[9.0, 2.5],
320                               end_point1=[13.0, 2.5],
321                               height=0.75,
322                               verbose=False)
323
324        domain.forcing_terms.append(culvert)
325
326
327        #-----------------------------------------------------------------------
328        # Setup boundary conditions
329        #-----------------------------------------------------------------------
330
331        # Inflow based on Flow Depth and Approaching Momentum
332
333        Br = Reflective_boundary(domain)              # Solid reflective wall
334        domain.set_boundary({'left': Br, 'right': Br, 'top': Br, 'bottom': Br})
335
336
337        #-----------------------------------------------------------------------
338        # Evolve system through time
339        #-----------------------------------------------------------------------
340
341        ref_volume = domain.get_quantity('stage').get_integral()
342        for t in domain.evolve(yieldstep = 1, finaltime = 25):
343            new_volume = domain.get_quantity('stage').get_integral()
344
345            msg = 'Total volume has changed'
346            assert num.allclose(new_volume, ref_volume, rtol=1.0e-10), msg
347            pass
348
349
350
351        os.remove('Test_culvert_dry.sww')
352
353    def test_that_culvert_flows_conserves_volume(self):
354        """test_that_culvert_flows_conserves_volume
355
356        Test that culvert on a sloping dry bed limits flows when very little water
357        is present at inlet.
358
359        Uses helper function: run_culvert_flow_problem(depth):
360
361        """
362
363        # Try this for a range of depths
364        for depth in [0.1, 0.2, 0.5, 1.0]:
365            run_culvert_flow_problem(depth)
366
367
368    def OBSOLETE_XXXtest_that_culvert_rating_limits_flow_in_shallow_inlet_condition(self):
369        """test_that_culvert_rating_limits_flow_in_shallow_inlet_condition
370
371        Test that culvert on a sloping dry bed limits flows when very little water
372        is present at inlet
373
374        This one is using the rating curve variant
375        """
376
377
378
379        path = get_pathname_from_package('anuga.culvert_flows')
380
381        length = 40.
382        width = 5.
383
384        dx = dy = 1           # Resolution: Length of subdivisions on both axes
385
386        points, vertices, boundary = rectangular_cross(int(length/dx),
387                                                       int(width/dy),
388                                                       len1=length,
389                                                       len2=width)
390        domain = Domain(points, vertices, boundary)
391        domain.set_name('Test_culvert_shallow') # Output name
392        domain.set_default_order(2)
393
394
395        #----------------------------------------------------------------------
396        # Setup initial conditions
397        #----------------------------------------------------------------------
398
399        def topography(x, y):
400            """Set up a weir
401
402            A culvert will connect either side
403            """
404            # General Slope of Topography
405            z=-x/1000
406
407            N = len(x)
408            for i in range(N):
409
410               # Sloping Embankment Across Channel
411                if 5.0 < x[i] < 10.1:
412                    # Cut Out Segment for Culvert face
413                    if  1.0+(x[i]-5.0)/5.0 <  y[i]  < 4.0 - (x[i]-5.0)/5.0:
414                       z[i]=z[i]
415                    else:
416                       z[i] +=  0.5*(x[i] -5.0)    # Sloping Segment  U/S Face
417                if 10.0 < x[i] < 12.1:
418                   z[i] +=  2.5                    # Flat Crest of Embankment
419                if 12.0 < x[i] < 14.5:
420                    # Cut Out Segment for Culvert face
421                    if  2.0-(x[i]-12.0)/2.5 <  y[i]  < 3.0 + (x[i]-12.0)/2.5:
422                       z[i]=z[i]
423                    else:
424                       z[i] +=  2.5-1.0*(x[i] -12.0) # Sloping D/S Face
425
426
427            return z
428
429
430        domain.set_quantity('elevation', topography)
431        domain.set_quantity('friction', 0.01)         # Constant friction
432        domain.set_quantity('stage',
433                            expression='elevation + 0.1') # Shallow initial condition
434
435        # Boyd culvert
436        culvert = Culvert_flow(domain,
437                               label='Culvert No. 1',
438                               description='This culvert is a test unit 1.2m Wide by 0.75m High',
439                               end_point0=[9.0, 2.5],
440                               end_point1=[13.0, 2.5],
441                               width=1.20, height=0.75,
442                               culvert_routine=boyd_generalised_culvert_model,
443                               number_of_barrels=1,
444                               update_interval=2,
445                               verbose=False)
446
447        # Rating curve
448        #filename = os.path.join(path, 'example_rating_curve.csv')
449        #culvert = Culvert_flow(domain,
450        #                       culvert_description_filename=filename,
451        #                       end_point0=[9.0, 2.5],
452        #                       end_point1=[13.0, 2.5],
453        #                       trigger_depth=0.01,
454        #                       verbose=False)
455
456        domain.forcing_terms.append(culvert)
457
458
459        #-----------------------------------------------------------------------
460        # Setup boundary conditions
461        #-----------------------------------------------------------------------
462
463        # Inflow based on Flow Depth and Approaching Momentum
464
465        Br = Reflective_boundary(domain)              # Solid reflective wall
466        domain.set_boundary({'left': Br, 'right': Br, 'top': Br, 'bottom': Br})
467
468
469
470        #-----------------------------------------------------------------------
471        # Evolve system through time
472        #-----------------------------------------------------------------------
473
474        print 'depth', 0.1
475        ref_volume = domain.get_quantity('stage').get_integral()
476        for t in domain.evolve(yieldstep = 0.1, finaltime = 25):
477            new_volume = domain.get_quantity('stage').get_integral()
478
479            msg = ('Total volume has changed: Is %.8f m^3 should have been %.8f m^3'
480                   % (new_volume, ref_volume))
481            assert num.allclose(new_volume, ref_volume, rtol=1.0e-10), msg
482
483
484        return
485        # Now try this again for a depth of 10 cm and for a range of other depths
486        for depth in [0.1, 0.2, 0.5, 1.0]:
487            print 'depth', depth
488            domain.set_time(0.0)
489
490            domain.set_quantity('elevation', topography)
491            domain.set_quantity('friction', 0.01)         # Constant friction
492            domain.set_quantity('stage',
493                                expression='elevation + %f' % depth)
494
495
496            ref_volume = domain.get_quantity('stage').get_integral()
497            for t in domain.evolve(yieldstep = 0.1, finaltime = 25):
498                new_volume = domain.get_quantity('stage').get_integral()
499
500                msg = 'Total volume has changed: Is %.8f m^3 should have been %.8f m^3'\
501                    % (new_volume, ref_volume)
502
503                assert num.allclose(new_volume, ref_volume, rtol=1.0e-10), msg
504
505
506
507    def test_that_culvert_dry_bed_boyd_does_not_produce_flow(self):
508        """test_that_culvert_in_dry_bed_boyd_does_not_produce_flow(self):
509
510        Test that culvert on a sloping dry bed doesn't produce flows
511        although there will be a 'pressure' head due to delta_w > 0
512
513        This one is using the 'Boyd' variant
514        """
515
516        path = get_pathname_from_package('anuga.culvert_flows')
517
518        length = 40.
519        width = 5.
520
521        dx = dy = 1           # Resolution: Length of subdivisions on both axes
522
523        points, vertices, boundary = rectangular_cross(int(length/dx),
524                                                       int(width/dy),
525                                                       len1=length,
526                                                       len2=width)
527        domain = Domain(points, vertices, boundary)
528        domain.set_name('Test_culvert_dry') # Output name
529        domain.set_default_order(2)
530
531
532        #----------------------------------------------------------------------
533        # Setup initial conditions
534        #----------------------------------------------------------------------
535
536        def topography(x, y):
537            """Set up a weir
538
539            A culvert will connect either side
540            """
541            # General Slope of Topography
542            z=-x/1000
543
544            N = len(x)
545            for i in range(N):
546
547               # Sloping Embankment Across Channel
548                if 5.0 < x[i] < 10.1:
549                    # Cut Out Segment for Culvert face
550                    if  1.0+(x[i]-5.0)/5.0 <  y[i]  < 4.0 - (x[i]-5.0)/5.0:
551                       z[i]=z[i]
552                    else:
553                       z[i] +=  0.5*(x[i] -5.0)    # Sloping Segment  U/S Face
554                if 10.0 < x[i] < 12.1:
555                   z[i] +=  2.5                    # Flat Crest of Embankment
556                if 12.0 < x[i] < 14.5:
557                    # Cut Out Segment for Culvert face
558                    if  2.0-(x[i]-12.0)/2.5 <  y[i]  < 3.0 + (x[i]-12.0)/2.5:
559                       z[i]=z[i]
560                    else:
561                       z[i] +=  2.5-1.0*(x[i] -12.0) # Sloping D/S Face
562
563
564            return z
565
566
567        domain.set_quantity('elevation', topography)
568        domain.set_quantity('friction', 0.01)         # Constant friction
569        domain.set_quantity('stage',
570                            expression='elevation')   # Dry initial condition
571
572
573        filename = os.path.join(path, 'example_rating_curve.csv')
574
575
576        culvert = Culvert_flow(domain,
577                               label='Culvert No. 1',
578                               description='This culvert is a test unit 1.2m Wide by 0.75m High',
579                               end_point0=[9.0, 2.5],
580                               end_point1=[13.0, 2.5],
581                               width=1.20, height=0.75,
582                               culvert_routine=boyd_generalised_culvert_model,
583                               number_of_barrels=1,
584                               update_interval=2,
585                               verbose=False)
586
587        domain.forcing_terms.append(culvert)
588
589
590        #-----------------------------------------------------------------------
591        # Setup boundary conditions
592        #-----------------------------------------------------------------------
593
594        # Inflow based on Flow Depth and Approaching Momentum
595
596        Br = Reflective_boundary(domain)              # Solid reflective wall
597        domain.set_boundary({'left': Br, 'right': Br, 'top': Br, 'bottom': Br})
598
599
600        #-----------------------------------------------------------------------
601        # Evolve system through time
602        #-----------------------------------------------------------------------
603
604        ref_volume = domain.get_quantity('stage').get_integral()
605        for t in domain.evolve(yieldstep = 1, finaltime = 25):
606
607            new_volume = domain.get_quantity('stage').get_integral()
608
609            msg = 'Total volume has changed'
610            assert num.allclose(new_volume, ref_volume, rtol=1.0e-10), msg
611            pass
612
613
614
615
616
617    def test_predicted_boyd_flow(self):
618        """test_predicted_boyd_flow
619
620        Test that flows predicted by the boyd method are consistent with what what
621        is calculated in engineering codes.
622        The data was supplied by Petar Milevski
623        """
624
625        # FIXME(Ole) this is nowhere near finished
626        path = get_pathname_from_package('anuga.culvert_flows')
627
628        length = 12.
629        width = 5.
630
631        dx = dy = 0.5           # Resolution: Length of subdivisions on both axes
632
633        points, vertices, boundary = rectangular_cross(int(length/dx),
634                                                       int(width/dy),
635                                                       len1=length,
636                                                       len2=width)
637        domain = Domain(points, vertices, boundary)
638
639        domain.set_name('test_culvert')                 # Output name
640        domain.set_default_order(2)
641
642
643        #----------------------------------------------------------------------
644        # Setup initial conditions
645        #----------------------------------------------------------------------
646
647        def topography(x, y):
648            # General Slope of Topography
649            z=-x/10
650
651            return z
652
653
654        domain.set_quantity('elevation', topography)
655        domain.set_quantity('friction', 0.01)         # Constant friction
656        domain.set_quantity('stage', expression='elevation')
657
658
659        Q0 = domain.get_quantity('stage')
660        Q1 = Quantity(domain)
661
662        # Add depths to stage
663        head_water_depth = 0.169
664        tail_water_depth = 0.089
665
666        inlet_poly = [[0,0], [6,0], [6,5], [0,5]]
667        outlet_poly = [[6,0], [12,0], [12,5], [6,5]]
668
670                                        (outlet_poly, tail_water_depth)]))
671
672        domain.set_quantity('stage', Q0 + Q1)
673
674
675
676        culvert = Culvert_flow(domain,
677                               label='Test culvert',
678                               description='4 m test culvert',
679                               end_point0=[4.0, 2.5],
680                               end_point1=[8.0, 2.5],
681                               width=1.20,
682                               height=0.75,
683                               culvert_routine=boyd_generalised_culvert_model,
684                               number_of_barrels=1,
685                               verbose=False)
686
687
688        domain.forcing_terms.append(culvert)
689
690        # Call
691        culvert(domain)
692
693
694
695
696    def test_momentum_jet(self):
697        """test_momentum_jet
698
699        Test that culvert_class can accept keyword use_momentum_jet
700        This does not yet imply that the values have been tested. FIXME
701        """
702
703
704        length = 40.
705        width = 5.
706
707        dx = dy = 1           # Resolution: Length of subdivisions on both axes
708
709        points, vertices, boundary = rectangular_cross(int(length/dx),
710                                                       int(width/dy),
711                                                       len1=length,
712                                                       len2=width)
713        domain = Domain(points, vertices, boundary)
714        domain.set_name('Test_culvert_shallow') # Output name
715        domain.set_default_order(2)
716
717
718        #----------------------------------------------------------------------
719        # Setup initial conditions
720        #----------------------------------------------------------------------
721
722        def topography(x, y):
723            """Set up a weir
724
725            A culvert will connect either side
726            """
727            # General Slope of Topography
728            z=-x/1000
729
730            N = len(x)
731            for i in range(N):
732
733               # Sloping Embankment Across Channel
734                if 5.0 < x[i] < 10.1:
735                    # Cut Out Segment for Culvert face
736                    if  1.0+(x[i]-5.0)/5.0 <  y[i]  < 4.0 - (x[i]-5.0)/5.0:
737                       z[i]=z[i]
738                    else:
739                       z[i] +=  0.5*(x[i] -5.0)    # Sloping Segment  U/S Face
740                if 10.0 < x[i] < 12.1:
741                   z[i] +=  2.5                    # Flat Crest of Embankment
742                if 12.0 < x[i] < 14.5:
743                    # Cut Out Segment for Culvert face
744                    if  2.0-(x[i]-12.0)/2.5 <  y[i]  < 3.0 + (x[i]-12.0)/2.5:
745                       z[i]=z[i]
746                    else:
747                       z[i] +=  2.5-1.0*(x[i] -12.0) # Sloping D/S Face
748
749
750            return z
751
752
753        domain.set_quantity('elevation', topography)
754        domain.set_quantity('friction', 0.01)         # Constant friction
755        domain.set_quantity('stage',
756                            expression='elevation + 1.0') # Shallow initial condition
757
758        # Boyd culvert
759        culvert = Culvert_flow(domain,
760                               label='Culvert No. 1',
761                               description='This culvert is a test unit 1.2m Wide by 0.75m High',
762                               end_point0=[9.0, 2.5],
763                               end_point1=[13.0, 2.5],
764                               width=1.20, height=0.75,
765                               culvert_routine=boyd_generalised_culvert_model,
766                               number_of_barrels=1,
767                               use_momentum_jet=True,
768                               update_interval=2,
769                               verbose=False)
770
771
772        domain.forcing_terms.append(culvert)
773
774
775        # Call
776        culvert(domain)
777
778
779        #-----------------------------------------------------------------------
780        # Setup boundary conditions
781        #-----------------------------------------------------------------------
782
783
784        Br = Reflective_boundary(domain)              # Solid reflective wall
785        domain.set_boundary({'left': Br, 'right': Br, 'top': Br, 'bottom': Br})
786
787        #-----------------------------------------------------------------------
788        # Evolve system through time
789        #-----------------------------------------------------------------------
790
791        ref_volume = domain.get_quantity('stage').get_integral()
792        for t in domain.evolve(yieldstep = 0.1, finaltime = 25):
793            new_volume = domain.get_quantity('stage').get_integral()
794
795            msg = ('Total volume has changed: Is %.8f m^3 should have been %.8f m^3'
796                   % (new_volume, ref_volume))
797            assert num.allclose(new_volume, ref_volume, rtol=1.0e-10), msg
798
799
800
801
802
803#-------------------------------------------------------------
804
805if __name__ == "__main__":
806    suite = unittest.makeSuite(Test_Culvert, 'test')
807    runner = unittest.TextTestRunner() #verbosity=2)
808    runner.run(suite)
809
Note: See TracBrowser for help on using the repository browser.