Changeset 4697
 Timestamp:
 Sep 3, 2007, 6:29:38 PM (16 years ago)
 File:

 1 edited
Legend:
 Unmodified
 Added
 Removed

anuga_work/development/Rudy_vandrie_2007/culvert_flow.py
r4440 r4697 12 12 # Import necessary modules 13 13 # 14 from anuga.abstract_2d_finite_volumes.mesh_factory import rectangular 14 from anuga.abstract_2d_finite_volumes.mesh_factory import rectangular_cross 15 15 from anuga.shallow_water import Domain 16 16 from anuga.shallow_water.shallow_water_domain import Reflective_boundary 17 17 from anuga.shallow_water.shallow_water_domain import Dirichlet_boundary 18 18 19 19 20 # 20 21 # Setup computational domain 21 22 # 22 23 #N = 120 24 N = 40 25 points, vertices, boundary = rectangular(N, N) 23 length = 40. 24 width = 5. 25 26 #dx = dy = 1 # Resolution: Length of subdivisions on both axes 27 #dx = dy = .5 # Resolution: Length of subdivisions on both axes 28 dx = dy = .1 # Resolution: Length of subdivisions on both axes 29 30 points, vertices, boundary = rectangular_cross(int(length/dx), int(width/dy), 31 len1=length, len2=width) 26 32 domain = Domain(points, vertices, boundary) 27 domain.set_name(' kitchen_sink') # Output name33 domain.set_name('culvert_example') # Output name 28 34 print 'Size', len(domain) 29 35 30 31 36 # 32 37 # Setup initial conditions 33 38 # 34 39 35 domain.set_quantity('elevation', 0.0) # Zero bed elevation 36 domain.set_quantity('stage', 0.015) # Constant stage 37 domain.set_quantity('friction', 0.005) # Constant friction 40 def topography(x, y): 41 """Set up a weir 42 43 A culvert will connect either side 44 """ 45 46 z = x/10 47 48 N = len(x) 49 for i in range(N): 50 51 # Weir from 10 to 12 m 52 if 10 < x[i] < 12: 53 z[i] += 0.4  0.05*y[i] 54 55 # Constriction 56 #if 27 < x[i] < 29 and y[i] > 3: 57 # z[i] += 2 58 59 # Pole 60 #if (x[i]  34)**2 + (y[i]  2)**2 < 0.4**2: 61 # z[i] += 2 62 63 return z 64 65 66 domain.set_quantity('elevation', topography) # Use function for elevation 67 domain.set_quantity('friction', 0.01) # Constant friction 68 domain.set_quantity('stage', 69 expression='elevation') # Dry initial condition 70 71 38 72 39 73 … … 67 101 # specified area 68 102 """ 69 103 70 104 # FIXME (OLE): Add a polygon as an alternative. 71 105 # FIXME (OLE): Generalise to all quantities 72 73 def __init__(self, 106 # FIXME (OLE): 107 108 def __init__(self, 109 domain, 74 110 center=None, radius=None, 75 111 flow=0.0, 76 112 quantity_name = 'stage'): 77 113 114 from anuga.utilities.numerical_tools import ensure_numeric 115 78 116 if center is not None and radius is not None: 79 117 assert len(center) == 2 … … 81 119 msg = 'Both center and radius must be specified' 82 120 raise Exception, msg 83 84 self.center = center 121 122 self.domain = domain 123 self.center = ensure_numeric(center) 85 124 self.radius = radius 86 125 self.flow = flow 126 87 127 self.quantity = domain.quantities[quantity_name].explicit_update 88 128 129 # Determine indices in flow area 130 N = len(domain) 131 self.indices = [] 132 coordinates = domain.get_centroid_coordinates() 133 for k in range(N): 134 x, y = coordinates[k,:] # Centroid 135 if ((xcenter[0])**2+(ycenter[1])**2) < radius**2: 136 self.indices.append(k) 137 89 138 90 139 def __call__(self, domain): 91 140 92 # Determine indices in flow area93 if not hasattr(self, 'indices'):94 center = self.center95 radius = self.radius96 97 N = len(domain)98 self.indices = []99 coordinates = domain.get_centroid_coordinates()100 for k in range(N):101 x, y = coordinates[k,:] # Centroid102 if ((xcenter[0])**2+(ycenter[1])**2) < radius**2:103 self.indices.append(k)104 105 141 # Update inflow 106 142 if callable(self.flow): … … 117 153 """ 118 154 119 def __init__(self, 120 inflow_center=None, inflow_radius=None, 121 outflow_center=None, outflow_radius=None, 122 transfer_flow=None) 123 124 self.inflow=Inflow(center=inflow_center, 125 radius=inflow_radius, 126 flow=transfer_flow) 155 def __init__(self, 156 domain, 157 center0=None, radius0=None, 158 center1=None, radius1=None, 159 transfer_flow=None): 160 161 from Numeric import sqrt, sum 162 163 self.openings = [] 164 self.openings.append(Inflow(domain, 165 center=center0, 166 radius=radius0, 167 flow=None)) 127 168 128 self.outflow=Inflow(center=outflow_center, 129 radius=outflow_radius, 130 flow=transfer_flow) #!!! 169 self.openings.append(Inflow(domain, 170 center=center1, 171 radius=radius1, 172 flow=None)) 131 173 132 174 # Compute distance between opening centers. I assume there are only two. 175 self.distance = sqrt(sum( (self.openings[0].center  self.openings[1].center)**2 )) 176 print 'Distance between openings is', self.distance 177 133 178 134 179 def __call__(self, domain): 135 self.inflow(domain) 136 self.outflow(domain) 180 from anuga.utilities.numerical_tools import mean 181 182 183 # Get average water depths at each opening 184 for opening in self.openings: 185 stage = domain.quantities['stage'].get_values(location='centroids', 186 indices=opening.indices) 187 elevation = domain.quantities['elevation'].get_values(location='centroids', 188 indices=opening.indices) 189 190 # Store current avg depth with opening object 191 opening.depth = mean(stageelevation) 192 # print 'Depth', opening.depth 193 194 195 # Here's a simple transfer function hardwired into this Culvert class: 196 # Let flow rate depend on difference in depth 197 # Flow is instantaneous for now but one might use 198 # the precomputed distance between the two openings 199 # to somehow introduce a delay. 200 201 delta_h = self.openings[0].depth  self.openings[1].depth 202 flow_rate = 0.5*9.81*delta_h**2 # Just an example of flow rate 203 flow_rate /= (self.openings[0].radius + self.openings[1].radius)/2 # 204 205 if delta_h > 0: 206 # Opening 0 has the greatest depth. 207 # Flow will go from opening 0 to opening 1 208 209 self.openings[0].flow = flow_rate 210 self.openings[1].flow = flow_rate 211 else: 212 # Opening 1 has the greatest depth. 213 # Flow will go from opening 1 to opening 0 214 215 self.openings[0].flow = flow_rate 216 self.openings[1].flow = flow_rate 217 218 219 # Execute flow term for each opening 220 for opening in self.openings: 221 opening(domain) 222 137 223 138 224 139 sink = Inflow(center=[0.7, 0.4], radius=0.0707, flow = 0.0) 140 tap = Inflow(center=(0.5, 0.5), radius=0.0316, 141 flow=lambda t: min(4*t, 5)) # Tap turning up 142 143 144 domain.forcing_terms.append(tap) 145 domain.forcing_terms.append(sink) 225 #sink = Inflow(domain, center=[0.7, 0.4], radius=0.0707, flow = 0.0) 226 #tap = Inflow(domain, center=(0.5, 0.5), radius=0.0316, 227 # flow=lambda t: min(4*t, 5)) # Tap turning up 228 #domain.forcing_terms.append(tap) 229 #domain.forcing_terms.append(sink) 230 231 culvert = Culvert_flow(domain, 232 center0=[5.0, 2.5], radius0=0.3, 233 center1=[20., 2.5], radius1=0.3, 234 transfer_flow=None) # Hardwired for now 235 236 domain.forcing_terms.append(culvert) 146 237 147 238 # 148 239 # Setup boundary conditions 149 240 # 150 241 Bi = Dirichlet_boundary([0.4, 0, 0]) # Inflow 151 242 Br = Reflective_boundary(domain) # Solid reflective wall 152 domain.set_boundary({'left': Br, 'right': Br, 'top': Br, 'bottom': Br}) 243 Bo = Dirichlet_boundary([5, 0, 0]) # Outflow 244 245 domain.set_boundary({'left': Bi, 'right': Bo, 'top': Br, 'bottom': Br}) 246 153 247 154 248 # … … 158 252 domain.write_time() 159 253 160 if domain.get_time() >= 4 and tap.flow != 0.0:161 print 'Turning tap off'162 tap.flow = 0.0254 #if domain.get_time() >= 4 and tap.flow != 0.0: 255 # print 'Turning tap off' 256 # tap.flow = 0.0 163 257 164 if domain.get_time() >= 3 and sink.flow < 0.0:165 print 'Turning drain on'166 sink.flow = 0.8167 258 #if domain.get_time() >= 3 and sink.flow < 0.0: 259 # print 'Turning drain on' 260 # sink.flow = 0.8 261
Note: See TracChangeset
for help on using the changeset viewer.