Changeset 6079


Ignore:
Timestamp:
Dec 15, 2008, 7:15:59 PM (15 years ago)
Author:
ole
Message:

Fixed dry bed culvert behaviour

Location:
anuga_core/source/anuga/culvert_flows
Files:
2 edited

Legend:

Unmodified
Added
Removed
  • anuga_core/source/anuga/culvert_flows/culvert_class.py

    r6059 r6079  
    308308                   
    309309                # Store current average stage and depth with each opening object
     310                opening.depth = stage - opening.elevation
    310311                opening.stage = stage
    311312
     
    327328
    328329
    329             # Calculate discharge for one barrel and set inlet.rate and outlet.rate
    330             try:
    331                 Q = interpolate_linearly(delta_w, self.rating_curve[:,0], self.rating_curve[:,1])
    332             except Below_interval, e:
    333                 Q = self.rating_curve[0,1]             
    334                 msg = '%.2fs: Delta head smaller than rating curve minimum: ' %time
    335                 msg += str(e)
    336                 msg += '\n        I will use minimum discharge %.2f m^3/s for culvert "%s"'\
    337                     %(Q, self.label)
    338                 if hasattr(self, 'log_filename'):                   
    339                     log_to_file(self.log_filename, msg)
    340             except Above_interval, e:
    341                 Q = self.rating_curve[-1,1]             
    342                 msg = '%.2fs: Delta head greater than rating curve maximum: ' %time
    343                 msg += str(e)
    344                 msg += '\n        I will use maximum discharge %.2f m^3/s for culvert "%s"'\
    345                     %(Q, self.label)
    346                 if hasattr(self, 'log_filename'):                   
    347                     log_to_file(self.log_filename, msg)                   
     330            if self.inlet.depth <= 0.01:
     331                Q = 0.0
     332            else:
     333                # Calculate discharge for one barrel and set inlet.rate and outlet.rate
     334                try:
     335                    Q = interpolate_linearly(delta_w, self.rating_curve[:,0], self.rating_curve[:,1])
     336                except Below_interval, e:
     337                    Q = self.rating_curve[0,1]             
     338                    msg = '%.2fs: Delta head smaller than rating curve minimum: ' %time
     339                    msg += str(e)
     340                    msg += '\n        I will use minimum discharge %.2f m^3/s for culvert "%s"'\
     341                        %(Q, self.label)
     342                    if hasattr(self, 'log_filename'):                   
     343                        log_to_file(self.log_filename, msg)
     344                except Above_interval, e:
     345                    Q = self.rating_curve[-1,1]             
     346                    msg = '%.2fs: Delta head greater than rating curve maximum: ' %time
     347                    msg += str(e)
     348                    msg += '\n        I will use maximum discharge %.2f m^3/s for culvert "%s"'\
     349                        %(Q, self.label)
     350                    if hasattr(self, 'log_filename'):                   
     351                        log_to_file(self.log_filename, msg)                   
    348352
    349353               
  • anuga_core/source/anuga/culvert_flows/test_culvert_class.py

    r6076 r6079  
    1212    Transmissive_boundary, Time_boundary
    1313
    14 from anuga.culvert_flows.culvert_class import Culvert_flow_rating
     14from anuga.culvert_flows.culvert_class import Culvert_flow_rating, Culvert_flow_energy
    1515from anuga.culvert_flows.culvert_routines import boyd_generalised_culvert_model
    1616     
     
    124124
    125125
    126     def test_that_culvert_dry_bed(self):
     126    def test_that_culvert_dry_bed_rating(self):
    127127        """test_that_culvert_in_dry_bed_does_not_produce_flow(self):
    128128       
    129129        Test that culvert on a sloping dry bed doesn't produce flows
    130130        although there will be a 'pressure' head due to delta_w > 0
     131
     132        This one is using the rating curve variant
    131133        """
    132134
     
    214216        ref_volume = domain.get_quantity('stage').get_integral()
    215217        for t in domain.evolve(yieldstep = 1, finaltime = 25):
    216             print domain.timestepping_statistics()
     218            #print domain.timestepping_statistics()
     219            new_volume = domain.get_quantity('stage').get_integral()
     220           
     221            msg = 'Total volume has changed'
     222            assert allclose(new_volume, ref_volume), msg
     223            pass
     224   
     225
     226
     227    def test_that_culvert_dry_bed_boyd(self):
     228        """test_that_culvert_in_dry_bed_does_not_produce_flow(self):
     229       
     230        Test that culvert on a sloping dry bed doesn't produce flows
     231        although there will be a 'pressure' head due to delta_w > 0
     232
     233        This one is using the 'Boyd' variant       
     234        """
     235
     236        path = get_pathname_from_package('anuga.culvert_flows')   
     237       
     238        length = 40.
     239        width = 5.
     240
     241        dx = dy = 1           # Resolution: Length of subdivisions on both axes
     242
     243        points, vertices, boundary = rectangular_cross(int(length/dx),
     244                                                       int(width/dy),
     245                                                       len1=length,
     246                                                       len2=width)
     247        domain = Domain(points, vertices, boundary)   
     248        domain.set_name('Test_culvert_dry') # Output name
     249        domain.set_default_order(2)
     250
     251
     252        #----------------------------------------------------------------------
     253        # Setup initial conditions
     254        #----------------------------------------------------------------------
     255
     256        def topography(x, y):
     257            """Set up a weir
     258           
     259            A culvert will connect either side
     260            """
     261            # General Slope of Topography
     262            z=-x/1000
     263           
     264            N = len(x)
     265            for i in range(N):
     266
     267               # Sloping Embankment Across Channel
     268                if 5.0 < x[i] < 10.1:
     269                    # Cut Out Segment for Culvert FACE               
     270                    if  1.0+(x[i]-5.0)/5.0 <  y[i]  < 4.0 - (x[i]-5.0)/5.0:
     271                       z[i]=z[i]
     272                    else:
     273                       z[i] +=  0.5*(x[i] -5.0)    # Sloping Segment  U/S Face
     274                if 10.0 < x[i] < 12.1:
     275                   z[i] +=  2.5                    # Flat Crest of Embankment
     276                if 12.0 < x[i] < 14.5:
     277                    # Cut Out Segment for Culvert FACE               
     278                    if  2.0-(x[i]-12.0)/2.5 <  y[i]  < 3.0 + (x[i]-12.0)/2.5:
     279                       z[i]=z[i]
     280                    else:
     281                       z[i] +=  2.5-1.0*(x[i] -12.0) # Sloping D/S Face
     282                                   
     283                       
     284            return z
     285
     286
     287        domain.set_quantity('elevation', topography)
     288        domain.set_quantity('friction', 0.01)         # Constant friction
     289        domain.set_quantity('stage',
     290                            expression='elevation')   # Dry initial condition
     291
     292
     293        filename = os.path.join(path, 'example_rating_curve.csv')
     294
     295
     296        culvert = Culvert_flow_energy(domain,
     297                                      label='Culvert No. 1',
     298                                      description='This culvert is a test unit 1.2m Wide by 0.75m High',   
     299                                      end_point0=[9.0, 2.5],
     300                                      end_point1=[13.0, 2.5],
     301                                      width=1.20,height=0.75,
     302                                      culvert_routine=boyd_generalised_culvert_model,       
     303                                      number_of_barrels=1,
     304                                      update_interval=2,
     305                                      verbose=True)
     306       
     307        domain.forcing_terms.append(culvert)
     308       
     309
     310        #-----------------------------------------------------------------------
     311        # Setup boundary conditions
     312        #-----------------------------------------------------------------------
     313
     314        # Inflow based on Flow Depth and Approaching Momentum
     315
     316        Br = Reflective_boundary(domain)              # Solid reflective wall
     317        domain.set_boundary({'left': Br, 'right': Br, 'top': Br, 'bottom': Br})
     318
     319
     320        #-----------------------------------------------------------------------
     321        # Evolve system through time
     322        #-----------------------------------------------------------------------
     323
     324        ref_volume = domain.get_quantity('stage').get_integral()
     325        for t in domain.evolve(yieldstep = 1, finaltime = 25):
     326            #print domain.timestepping_statistics()
    217327            new_volume = domain.get_quantity('stage').get_integral()
    218328           
Note: See TracChangeset for help on using the changeset viewer.