[5773] | 1 | """ Testing CULVERT (Changing from Horizontal Abstraction to Vertical Abstraction |
---|
| 2 | |
---|
| 3 | This example includes a Model Topography that shows a TYPICAL Headwall Configuration |
---|
| 4 | |
---|
| 5 | The aim is to change the Culvert Routine to Model more precisely the abstraction |
---|
| 6 | from a vertical face. |
---|
| 7 | |
---|
| 8 | The inflow must include the impact of Approach velocity. |
---|
| 9 | Similarly the Outflow has MOMENTUM Not just Up welling as in the Horizontal Style |
---|
| 10 | abstraction |
---|
| 11 | |
---|
| 12 | """ |
---|
| 13 | print 'Starting.... Importing Modules...' |
---|
| 14 | #------------------------------------------------------------------------------ |
---|
| 15 | # Import necessary modules |
---|
| 16 | #------------------------------------------------------------------------------ |
---|
| 17 | from anuga.abstract_2d_finite_volumes.mesh_factory import rectangular_cross |
---|
| 18 | |
---|
| 19 | from anuga.shallow_water import Domain, Reflective_boundary,\ |
---|
| 20 | Dirichlet_boundary,\ |
---|
| 21 | Transmissive_boundary, Time_boundary |
---|
| 22 | |
---|
[6143] | 23 | from anuga.culvert_flows.culvert_class import Culvert_flow |
---|
[5774] | 24 | from anuga.culvert_flows.culvert_routines import boyd_generalised_culvert_model |
---|
| 25 | |
---|
[5773] | 26 | from math import pi,pow,sqrt |
---|
[6150] | 27 | |
---|
[6304] | 28 | import numpy as num |
---|
[6150] | 29 | |
---|
| 30 | |
---|
[5773] | 31 | #------------------------------------------------------------------------------ |
---|
| 32 | # Setup computational domain |
---|
| 33 | #------------------------------------------------------------------------------ |
---|
| 34 | print 'Setting up domain' |
---|
| 35 | |
---|
| 36 | length = 40. |
---|
| 37 | width = 5. |
---|
| 38 | |
---|
[5777] | 39 | dx = dy = 1 # Resolution: Length of subdivisions on both axes |
---|
[5773] | 40 | #dx = dy = .5 # Resolution: Length of subdivisions on both axes |
---|
[5777] | 41 | #dx = dy = .5 # Resolution: Length of subdivisions on both axes |
---|
[5773] | 42 | #dx = dy = .1 # Resolution: Length of subdivisions on both axes |
---|
| 43 | |
---|
| 44 | points, vertices, boundary = rectangular_cross(int(length/dx), int(width/dy), |
---|
| 45 | len1=length, len2=width) |
---|
| 46 | domain = Domain(points, vertices, boundary) |
---|
| 47 | domain.set_name('Test_Culv_Flat_WL') # Output name |
---|
| 48 | domain.set_default_order(2) |
---|
| 49 | domain.H0 = 0.01 |
---|
| 50 | domain.tight_slope_limiters = 1 |
---|
| 51 | |
---|
| 52 | print 'Size', len(domain) |
---|
| 53 | |
---|
| 54 | #------------------------------------------------------------------------------ |
---|
| 55 | # Setup initial conditions |
---|
| 56 | #------------------------------------------------------------------------------ |
---|
| 57 | |
---|
| 58 | def topography(x, y): |
---|
| 59 | """Set up a weir |
---|
| 60 | |
---|
| 61 | A culvert will connect either side |
---|
| 62 | """ |
---|
| 63 | # General Slope of Topography |
---|
| 64 | z=-x/1000 |
---|
| 65 | |
---|
[5774] | 66 | # NOW Add bits and Pieces to topography |
---|
[5773] | 67 | N = len(x) |
---|
| 68 | for i in range(N): |
---|
| 69 | |
---|
| 70 | # Sloping Embankment Across Channel |
---|
| 71 | if 5.0 < x[i] < 10.1: |
---|
| 72 | if 1.0+(x[i]-5.0)/5.0 < y[i] < 4.0 - (x[i]-5.0)/5.0: # Cut Out Segment for Culvert FACE |
---|
| 73 | z[i]=z[i] |
---|
| 74 | else: |
---|
| 75 | z[i] += 0.5*(x[i] -5.0) # Sloping Segment U/S Face |
---|
| 76 | if 10.0 < x[i] < 12.1: |
---|
| 77 | z[i] += 2.5 # Flat Crest of Embankment |
---|
| 78 | if 12.0 < x[i] < 14.5: |
---|
| 79 | if 2.0-(x[i]-12.0)/2.5 < y[i] < 3.0 + (x[i]-12.0)/2.5: # Cut Out Segment for Culvert FACE |
---|
| 80 | z[i]=z[i] |
---|
| 81 | else: |
---|
| 82 | z[i] += 2.5-1.0*(x[i] -12.0) # Sloping D/S Face |
---|
| 83 | |
---|
| 84 | |
---|
| 85 | |
---|
| 86 | return z |
---|
| 87 | |
---|
| 88 | print 'Setting Quantities....' |
---|
| 89 | domain.set_quantity('elevation', topography) # Use function for elevation |
---|
| 90 | domain.set_quantity('friction', 0.01) # Constant friction |
---|
| 91 | domain.set_quantity('stage', |
---|
| 92 | expression='elevation') # Dry initial condition |
---|
| 93 | |
---|
| 94 | |
---|
| 95 | |
---|
| 96 | |
---|
| 97 | #------------------------------------------------------------------------------ |
---|
| 98 | # Setup specialised forcing terms |
---|
| 99 | #------------------------------------------------------------------------------ |
---|
| 100 | |
---|
| 101 | #------------------------------------------------------------------------------ |
---|
| 102 | # Setup CULVERT INLETS and OUTLETS in Current Topography |
---|
| 103 | #------------------------------------------------------------------------------ |
---|
| 104 | print 'DEFINING any Structures if Required' |
---|
| 105 | |
---|
| 106 | # DEFINE CULVERT INLET AND OUTLETS |
---|
| 107 | |
---|
[5774] | 108 | |
---|
[6143] | 109 | culvert_rating = Culvert_flow(domain, |
---|
[6055] | 110 | culvert_description_filename='example_rating_curve.csv', |
---|
[5774] | 111 | end_point0=[9.0, 2.5], |
---|
| 112 | end_point1=[13.0, 2.5], |
---|
| 113 | verbose=True) |
---|
| 114 | |
---|
[6055] | 115 | |
---|
[6143] | 116 | culvert_energy = Culvert_flow(domain, |
---|
[6059] | 117 | label='Culvert No. 1', |
---|
| 118 | description='This culvert is a test unit 1.2m Wide by 0.75m High', |
---|
| 119 | end_point0=[9.0, 2.5], |
---|
| 120 | end_point1=[13.0, 2.5], |
---|
| 121 | width=1.20,height=0.75, |
---|
| 122 | culvert_routine=boyd_generalised_culvert_model, |
---|
| 123 | number_of_barrels=1, |
---|
| 124 | update_interval=2, |
---|
[6517] | 125 | log_file=True, |
---|
| 126 | discharge_hydrograph=True, |
---|
[6059] | 127 | verbose=True) |
---|
[6055] | 128 | |
---|
[6059] | 129 | domain.forcing_terms.append(culvert_energy) |
---|
[5774] | 130 | |
---|
[5773] | 131 | #------------------------------------------------------------------------------ |
---|
| 132 | # Setup boundary conditions |
---|
| 133 | #------------------------------------------------------------------------------ |
---|
| 134 | print 'Setting Boundary Conditions' |
---|
| 135 | Bi = Dirichlet_boundary([0.0, 0.0, 0.0]) # Inflow based on Flow Depth and Approaching Momentum !!! |
---|
| 136 | Br = Reflective_boundary(domain) # Solid reflective wall |
---|
| 137 | Bo = Dirichlet_boundary([-5, 0, 0]) # Outflow |
---|
[6150] | 138 | Btus = Time_boundary(domain, lambda t: [0.0+ 1.25*(1+num.sin(2*pi*(t-4)/10)), 0.0, 0.0]) |
---|
| 139 | Btds = Time_boundary(domain, lambda t: [0.0+ 0.75*(1+num.sin(2*pi*(t-4)/20)), 0.0, 0.0]) |
---|
[5773] | 140 | domain.set_boundary({'left': Btus, 'right': Btds, 'top': Br, 'bottom': Br}) |
---|
| 141 | |
---|
| 142 | |
---|
| 143 | #------------------------------------------------------------------------------ |
---|
| 144 | # Evolve system through time |
---|
| 145 | #------------------------------------------------------------------------------ |
---|
| 146 | |
---|
[6051] | 147 | #for t in domain.evolve(yieldstep = 1, finaltime = 25): |
---|
| 148 | # print domain.timestepping_statistics() |
---|
[5862] | 149 | |
---|
[5773] | 150 | |
---|
[5864] | 151 | |
---|
| 152 | |
---|
[5868] | 153 | #import sys; sys.exit() |
---|
[5862] | 154 | # Profiling code |
---|
| 155 | import time |
---|
| 156 | t0 = time.time() |
---|
[5773] | 157 | |
---|
[6051] | 158 | s = 'for t in domain.evolve(yieldstep = 1, finaltime = 25): domain.write_time()' |
---|
[5773] | 159 | |
---|
[5862] | 160 | import profile, pstats |
---|
| 161 | FN = 'profile.dat' |
---|
[5773] | 162 | |
---|
[5862] | 163 | profile.run(s, FN) |
---|
| 164 | |
---|
| 165 | print 'That took %.2f seconds' %(time.time()-t0) |
---|
| 166 | |
---|
| 167 | S = pstats.Stats(FN) |
---|
| 168 | #S.sort_stats('time').print_stats(20) |
---|
| 169 | s = S.sort_stats('cumulative').print_stats(30) |
---|
| 170 | |
---|
| 171 | print s |
---|