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

Last change on this file since 7497 was 7497, checked in by ole, 15 years ago

Tighter tolerance in allclose for volumetric check (Stephen Roberts)

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