| 1 | """Simple water flow example using ANUGA |
|---|
| 2 | Water flowing down a channel with a floodplain |
|---|
| 3 | """ |
|---|
| 4 | |
|---|
| 5 | #------------------------------------------------------------------------------ |
|---|
| 6 | # Import necessary modules |
|---|
| 7 | #------------------------------------------------------------------------------ |
|---|
| 8 | # Import standard shallow water domain and standard boundaries. |
|---|
| 9 | import anuga |
|---|
| 10 | import numpy |
|---|
| 11 | from anuga.structures.inlet_operator import Inlet_operator |
|---|
| 12 | from anuga import * |
|---|
| 13 | #from swb_domain import domain |
|---|
| 14 | #from anuga import * |
|---|
| 15 | #from balanced_basic import * |
|---|
| 16 | from balanced_dev import * |
|---|
| 17 | #from balanced_basic.swb2_domain import * |
|---|
| 18 | #from balanced_basic.swb2_domain import Domain |
|---|
| 19 | #------------------------------------------------------------------------------ |
|---|
| 20 | # Useful parameters for controlling this case |
|---|
| 21 | #------------------------------------------------------------------------------ |
|---|
| 22 | |
|---|
| 23 | floodplain_length = 1000.0 # Model domain length |
|---|
| 24 | floodplain_width = 14.0 # Model domain width |
|---|
| 25 | floodplain_slope = 1./300. |
|---|
| 26 | chan_initial_depth = 0.65 # Initial depth of water in the channel |
|---|
| 27 | chan_bankfull_depth = 1.0 # Bankfull depth of the channel |
|---|
| 28 | chan_width = 10.0 # Bankfull width of the channel |
|---|
| 29 | bankwidth = 2. # Width of the bank regions -- note that these protrude into the channel |
|---|
| 30 | man_n=0.03 # Manning's n |
|---|
| 31 | l0 = 1.000 # Length scale associated with triangle side length in channel (min_triangle area = 0.5*l0^2) |
|---|
| 32 | |
|---|
| 33 | assert chan_width < floodplain_width, \ |
|---|
| 34 | ' ERROR: Channel width is greater than floodplain width' |
|---|
| 35 | |
|---|
| 36 | assert bankwidth < chan_width/2., \ |
|---|
| 37 | 'ERROR: The bank width must be less than half the channel width' |
|---|
| 38 | |
|---|
| 39 | #------------------------------------------------------------------------------ |
|---|
| 40 | # Setup computational domain |
|---|
| 41 | #------------------------------------------------------------------------------ |
|---|
| 42 | |
|---|
| 43 | # Define boundary polygon -- in this case clockwise around the proposed boundary |
|---|
| 44 | boundary_polygon = [ [0.,0.], |
|---|
| 45 | [0., floodplain_length], |
|---|
| 46 | [floodplain_width/2. - chan_width/2., floodplain_length], |
|---|
| 47 | [floodplain_width/2. + chan_width/2., floodplain_length], |
|---|
| 48 | [floodplain_width, floodplain_length], |
|---|
| 49 | [floodplain_width, 0.], |
|---|
| 50 | [floodplain_width/2. + chan_width/2., 0.], |
|---|
| 51 | [floodplain_width/2. - chan_width/2., 0.] |
|---|
| 52 | ] |
|---|
| 53 | |
|---|
| 54 | # Define channel polygon, where we can optionally refine the resolution. |
|---|
| 55 | # Note that we set this a distance l0 inside the boundary polygon, so the polygons |
|---|
| 56 | # do not overlap. |
|---|
| 57 | channel_polygon = [ [floodplain_width/2. - chan_width/2., +l0], |
|---|
| 58 | [floodplain_width/2. - chan_width/2., floodplain_length-l0], |
|---|
| 59 | [floodplain_width/2. + chan_width/2., floodplain_length-l0], |
|---|
| 60 | [floodplain_width/2. + chan_width/2., +l0] |
|---|
| 61 | ] |
|---|
| 62 | |
|---|
| 63 | |
|---|
| 64 | # Define domain with appropriate boundary conditions |
|---|
| 65 | domain = create_domain_from_regions( boundary_polygon, |
|---|
| 66 | boundary_tags={'left': [0], |
|---|
| 67 | 'top1': [1], |
|---|
| 68 | 'chan_out': [2], |
|---|
| 69 | 'top2': [3], |
|---|
| 70 | 'right': [4], |
|---|
| 71 | 'bottom1': [5], |
|---|
| 72 | 'chan_in': [6], |
|---|
| 73 | 'bottom2': [7] }, |
|---|
| 74 | maximum_triangle_area = 0.5*l0*l0, |
|---|
| 75 | minimum_triangle_angle = 28.0, |
|---|
| 76 | mesh_filename = 'channel_floodplain1.msh', |
|---|
| 77 | interior_regions = [ ], |
|---|
| 78 | #interior_regions = [\ |
|---|
| 79 | # [channel_polygon, 0.5*l0*l0] ], |
|---|
| 80 | use_cache=False, |
|---|
| 81 | verbose=True) |
|---|
| 82 | |
|---|
| 83 | |
|---|
| 84 | domain.set_name('channel_floodplain1_bal_dev_lowvisc') # Output name |
|---|
| 85 | #domain.set_store_vertices_uniquely(True) |
|---|
| 86 | #domain.use_edge_limiter=False |
|---|
| 87 | #domain.extrapolate_velocity_second_order=False |
|---|
| 88 | #------------------------------------------------------------------------------ |
|---|
| 89 | # |
|---|
| 90 | # Setup initial conditions |
|---|
| 91 | # |
|---|
| 92 | #------------------------------------------------------------------------------ |
|---|
| 93 | |
|---|
| 94 | # Function for topography |
|---|
| 95 | def topography(x, y): |
|---|
| 96 | # Longitudinally sloping floodplain with channel in centre |
|---|
| 97 | #return -y*floodplain_slope -chan_bankfull_depth*\ |
|---|
| 98 | # (x>(floodplain_width/2. - chan_width/2.) )*\ |
|---|
| 99 | # (x<(floodplain_width/2. + chan_width/2.) ) |
|---|
| 100 | |
|---|
| 101 | elev1= -y*floodplain_slope - chan_bankfull_depth*\ |
|---|
| 102 | (x>(floodplain_width/2. - chan_width/2.))*\ |
|---|
| 103 | (x<(floodplain_width/2. + chan_width/2.)) |
|---|
| 104 | # Add banks |
|---|
| 105 | if(bankwidth>0.0): |
|---|
| 106 | leftbnk = floodplain_width/2. - chan_width/2. |
|---|
| 107 | rightbnk = floodplain_width/2. + chan_width/2. |
|---|
| 108 | # Left bank |
|---|
| 109 | elev2 = elev1 + (chan_bankfull_depth \ |
|---|
| 110 | - chan_bankfull_depth/bankwidth*(x - leftbnk))*\ |
|---|
| 111 | (x>leftbnk)*(x < leftbnk + bankwidth) |
|---|
| 112 | # Right bank |
|---|
| 113 | elev2 = elev2 + (chan_bankfull_depth \ |
|---|
| 114 | + chan_bankfull_depth/bankwidth*(x - rightbnk))*\ |
|---|
| 115 | (x>rightbnk-bankwidth)*(x < rightbnk) |
|---|
| 116 | |
|---|
| 117 | if(bankwidth==0.0): |
|---|
| 118 | elev2 = elev1 |
|---|
| 119 | # |
|---|
| 120 | return elev2 |
|---|
| 121 | |
|---|
| 122 | # |
|---|
| 123 | |
|---|
| 124 | #Function for stage |
|---|
| 125 | def stagetopo(x,y): |
|---|
| 126 | return -y*floodplain_slope -chan_bankfull_depth + chan_initial_depth |
|---|
| 127 | |
|---|
| 128 | domain.set_quantity('elevation', topography) # Use function for elevation |
|---|
| 129 | domain.set_quantity('friction', man_n) # Constant friction |
|---|
| 130 | domain.set_quantity('stage', stagetopo) # Use function for stage |
|---|
| 131 | domain.set_minimum_allowed_height(0.01) # |
|---|
| 132 | #def steady_uh_calc(): |
|---|
| 133 | # # Calculate a boundary condition potentially suitable for steady state flows. |
|---|
| 134 | # uh = (chan_initial_depth**(4./3.)*man_n**(-2.)*floodplain_slope)**0.5\ |
|---|
| 135 | # *chan_initial_depth # (u) * h |
|---|
| 136 | # return uh |
|---|
| 137 | # |
|---|
| 138 | #momin = steady_uh_calc() |
|---|
| 139 | #velin = momin/chan_initial_depth |
|---|
| 140 | # |
|---|
| 141 | #print 'momentum inflow is = ' + str(momin) |
|---|
| 142 | |
|---|
| 143 | # Define inlet operator |
|---|
| 144 | flow_in_yval=100.0 |
|---|
| 145 | if True: |
|---|
| 146 | line1 = [ [floodplain_width/2. - chan_width/2., flow_in_yval],\ |
|---|
| 147 | [floodplain_width/2. + chan_width/2., flow_in_yval] \ |
|---|
| 148 | ] |
|---|
| 149 | Qin = 0.5*(floodplain_slope*(chan_width*chan_initial_depth)**2.*man_n**(-2.)\ |
|---|
| 150 | *chan_initial_depth**(4./3.) )**0.5 |
|---|
| 151 | |
|---|
| 152 | Inlet_operator(domain, line1, Qin) |
|---|
| 153 | |
|---|
| 154 | print 'Discharge in = ', Qin #,'Velocity at inlet should be ', Qin/(chan_width*chan_initial_depth), \ |
|---|
| 155 | #'for rectangular channels, and approximately that for trapezoidal ones' |
|---|
| 156 | #------------------------------------------------------------------------------ |
|---|
| 157 | # |
|---|
| 158 | # Setup boundary conditions |
|---|
| 159 | # |
|---|
| 160 | #------------------------------------------------------------------------------ |
|---|
| 161 | |
|---|
| 162 | # We wish to define one inflow boundary at the upstream head of the channel, |
|---|
| 163 | # one at the outflow of the channel (e.g. transmissive) |
|---|
| 164 | # one at the floodplain edges, excluding the downstream floodplain edge |
|---|
| 165 | #(e.g. reflective) and one at the downstream floodplain edge (e.g. transmissive) |
|---|
| 166 | |
|---|
| 167 | #Bin_sub = anuga.shallow_water.boundaries.Dirichlet_discharge_boundary(domain,\ |
|---|
| 168 | # chan_initial_depth - chan_bankfull_depth , 0.) # An inflow for subcritical |
|---|
| 169 | #flow |
|---|
| 170 | #Bin_sup = anuga.shallow_water.boundaries.Dirichlet_discharge_boundary(domain, \ |
|---|
| 171 | # chan_initial_depth - chan_bankfull_depth , steady_uh_calc()) # An inflow |
|---|
| 172 | # # for |
|---|
| 173 | # #supercritical |
|---|
| 174 | # # flow |
|---|
| 175 | #Bin_sub = anuga.Dirichlet_boundary([chan_initial_depth - chan_bankfull_depth , \ |
|---|
| 176 | # 0., 0.]) # An inflow for subcritical flow -- maybe should use transmissive |
|---|
| 177 | # boundary instead |
|---|
| 178 | #Bin_sup = anuga.Dirichlet_boundary([chan_initial_depth - chan_bankfull_depth ,\ |
|---|
| 179 | # 0., steady_uh_calc()]) # An inflow for supercritical flow |
|---|
| 180 | |
|---|
| 181 | # Trial a novel boundary, which variously imposes 2 or 3 conserved variables |
|---|
| 182 | # depending on the local characteristic speed |
|---|
| 183 | #Bin_trial = anuga.shallow_water.boundaries.\ |
|---|
| 184 | # Dirichlet_inflow_boundary_sub_or_super_critical(domain, \ |
|---|
| 185 | # stage0=chan_initial_depth - chan_bankfull_depth, wh0 = momin,\ |
|---|
| 186 | # href = chan_initial_depth) |
|---|
| 187 | |
|---|
| 188 | # Trial a set_stage_transmissive_momentum inflow boundary |
|---|
| 189 | #def inflow_stage_boundary(t): |
|---|
| 190 | # return chan_initial_depth - chan_bankfull_depth |
|---|
| 191 | |
|---|
| 192 | # transmissive_momentum_set_stage |
|---|
| 193 | #Bin_tmss = anuga.shallow_water.boundaries.Transmissive_momentum_set_stage_boundary(domain=domain, function = inflow_stage_boundary) |
|---|
| 194 | |
|---|
| 195 | Br = anuga.Reflective_boundary(domain) # Solid reflective wall |
|---|
| 196 | Bt = anuga.Transmissive_boundary(domain) # Transmissive boundary |
|---|
| 197 | |
|---|
| 198 | |
|---|
| 199 | Bout_sub = anuga.Dirichlet_boundary( \ |
|---|
| 200 | [-floodplain_length*floodplain_slope - chan_bankfull_depth + \ |
|---|
| 201 | chan_initial_depth, 0., 0.]) #An outflow boundary for subcritical steady flow |
|---|
| 202 | |
|---|
| 203 | def outflow_stage_boundary(t): |
|---|
| 204 | return -floodplain_length*floodplain_slope \ |
|---|
| 205 | + chan_initial_depth - chan_bankfull_depth |
|---|
| 206 | |
|---|
| 207 | # Note that the outflow boundary may be slightly incorrect for the trapezoidal channel case, |
|---|
| 208 | # or incorrect more generally if there are numerical problems. But, in the central regions of |
|---|
| 209 | # the channel, this shouldn't really prevent us reaching steady, uniform flow. |
|---|
| 210 | Bout_tmss = anuga.shallow_water.boundaries.Transmissive_momentum_set_stage_boundary(domain, function = outflow_stage_boundary) |
|---|
| 211 | |
|---|
| 212 | domain.set_boundary({'left': Br, |
|---|
| 213 | 'right': Br, |
|---|
| 214 | 'top1': Bt, |
|---|
| 215 | 'top2': Bt, |
|---|
| 216 | 'bottom1': Br, |
|---|
| 217 | 'bottom2': Br, |
|---|
| 218 | 'chan_out': Bout_tmss, |
|---|
| 219 | 'chan_in': Br}) |
|---|
| 220 | |
|---|
| 221 | |
|---|
| 222 | # To monitor the boundary conditions, let's find centroid points representing sub and super critical flows. |
|---|
| 223 | #inpt = numpy.argmin( ( abs(domain.centroid_coordinates[:,0] - floodplain_width/2.) + \ |
|---|
| 224 | # abs(domain.centroid_coordinates[:,1] - 0.0) ) ) |
|---|
| 225 | #outpt = numpy.argmin( ( abs(domain.centroid_coordinates[:, 0] - floodplain_width/2.) +\ |
|---|
| 226 | # abs(domain.centroid_coordinates[:, 1] - floodplain_length) ) ) |
|---|
| 227 | |
|---|
| 228 | #------------------------------------------------------------------------------ |
|---|
| 229 | # |
|---|
| 230 | # Evolve system through time |
|---|
| 231 | # |
|---|
| 232 | #------------------------------------------------------------------------------ |
|---|
| 233 | |
|---|
| 234 | for t in domain.evolve(yieldstep=2.0, finaltime=3200.0): |
|---|
| 235 | print domain.timestepping_statistics() |
|---|
| 236 | xx=domain.quantities['ymomentum'].centroid_values |
|---|
| 237 | dd=(domain.quantities['stage'].centroid_values - domain.quantities['elevation'].centroid_values) |
|---|
| 238 | dd=dd*(dd>0.) |
|---|
| 239 | |
|---|
| 240 | tmp = xx/(dd+1.0e-06)*(dd>0.0) |
|---|
| 241 | print tmp.max(), tmp.argmax(), tmp.min(), tmp.argmin() |
|---|
| 242 | |
|---|
| 243 | # Compute flow through cross-section -- check that the inflow boundary condition is doing its job |
|---|
| 244 | # This also provides another useful steady-state check |
|---|
| 245 | if( numpy.floor(t/100.) == t/100. ): |
|---|
| 246 | print '#### COMPUTING FLOW THROUGH CROSS-SECTIONS########' |
|---|
| 247 | s1 = domain.get_flow_through_cross_section([[0., floodplain_length-300.0], [floodplain_width, floodplain_length-300.0]]) |
|---|
| 248 | s2 = domain.get_flow_through_cross_section([[0., floodplain_length-1.0], [floodplain_width, floodplain_length-1.0]]) |
|---|
| 249 | |
|---|
| 250 | print 'Cross sectional flows: ', s1, s2 |
|---|
| 251 | print '$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$' |
|---|
| 252 | |
|---|
| 253 | #vv = domain.get_flow_through_cross_section |
|---|
| 254 | #print domain.quantities['ymomentum'].get_integral() |
|---|
| 255 | #print 'Qin = ', Qin |
|---|
| 256 | #print domain.quantities['stage'].centroid_values[inpt] |
|---|
| 257 | #print domain.quantities['xmomentum'].centroid_values[inpt] |
|---|
| 258 | #print domain.quantities['ymomentum'].centroid_values[inpt] |
|---|
| 259 | |
|---|
| 260 | # Alter inflow boundary for sub and super critical flow |
|---|
| 261 | # FIXME: Varying the inflow boundary conditions does not work well with the modified ANUGA |
|---|
| 262 | # (but works great with the standard ANUGA). Especially the supercritical inflow seems bad. |
|---|
| 263 | # With the modified ANUGA, they oscillate (I wonder if the results are better with a finer grid?) |
|---|
| 264 | #h= domain.quantities['stage'].centroid_values[inpt] - domain.quantities['elevation'].centroid_values[inpt] |
|---|
| 265 | #uh=domain.quantities['xmomentum'].centroid_values[inpt] |
|---|
| 266 | #if( domain.g*h < (uh/h)**2): |
|---|
| 267 | # domain.set_boundary({'chan_in': Bin_sup}) |
|---|
| 268 | #else: |
|---|
| 269 | # domain.set_boundary({'chan_in': Bin_sub}) |
|---|
| 270 | |
|---|
| 271 | # Alter outflow boundary for sub and super critical flow |
|---|
| 272 | #h=domain.quantities['stage'].centroid_values[outpt] - domain.quantities['elevation'].centroid_values[outpt] |
|---|
| 273 | #uh=domain.quantities['xmomentum'].centroid_values[outpt] |
|---|
| 274 | #if( domain.g*h < (uh/h)**2): |
|---|
| 275 | # domain.set_boundary({'chan_out': Bt}) |
|---|
| 276 | #else: |
|---|
| 277 | # domain.set_boundary({'chan_out': Bout_sub}) |
|---|
| 278 | |
|---|
