Ignore:
Timestamp:
Jul 8, 2014, 12:59:46 PM (10 years ago)
Author:
davies
Message:

Code to read list of riverwall input files

File:
1 edited

Legend:

Unmodified
Added
Removed

  • trunk/anuga_work/development/gareth/tests/channel_floodplain/channel_floodplain1.py

    r9056 r9252  
    1010import numpy
    1111from anuga.structures.inlet_operator import Inlet_operator
    12 
    1312from anuga import create_domain_from_regions
    14 #from anuga import *
    15 #from swb_domain import domain
    16 #from anuga import *
    17 #from balanced_basic import *
    18 #from balanced_dev import *
    19 #from bal_and import *
    20 #from anuga_tsunami import *
    21 #from balanced_basic.swb2_domain import *
    22 #from balanced_basic.swb2_domain import Domain
    2313#------------------------------------------------------------------------------
    2414# Useful parameters for controlling this case
    … …  
    3020chan_initial_depth = 0.65 # Initial depth of water in the channel
    3121chan_bankfull_depth = 1.0 # Bankfull depth of the channel
    32 chan_width = 10.0 # Bankfull width of the channel
     22chan_width = 10. # Bankfull width of the channel
    3323bankwidth = 2. # Width of the bank regions -- note that these protrude into the channel
    3424man_n=0.03 # Manning's n
    35 l0 = 1.000 # Length scale associated with triangle side length in channel (min_triangle area = 0.5*l0^2)
     25l0 = 2.000 # Length scale associated with triangle side length in channel (min_triangle area = 0.5*l0^2)
    3626
    3727assert chan_width < floodplain_width, \
    … …  
    5646                     ]
    5747
    58 # Define channel polygon, where we can optionally refine the resolution.
    59 # Note that we set this a distance l0 inside the boundary polygon, so the polygons
    60 # do not overlap.
    61 channel_polygon = [ [floodplain_width/2. - chan_width/2., +l0],
    62                     [floodplain_width/2. - chan_width/2., floodplain_length-l0],
    63                     [floodplain_width/2. + chan_width/2., floodplain_length-l0],
    64                     [floodplain_width/2. + chan_width/2., +l0]
    65                     ]
     48## Define channel polygon, where we can optionally refine the resolution.
     49## Note that we set this a distance l0 inside the boundary polygon, so the polygons
     50## do not overlap.
    6651
     52breakLines={'n1':[[floodplain_width/2.- chan_width/2., floodplain_length],
     53                 [floodplain_width/2. - chan_width/2.,0.]],
     54            'n2':[[floodplain_width/2.+ chan_width/2., floodplain_length],
     55                  [floodplain_width/2. + chan_width/2.,0.]],
     56            # These ones are inside the channel
     57            'n3':[[floodplain_width/2.- chan_width/2.+bankwidth, floodplain_length],
     58                 [floodplain_width/2. - chan_width/2.+bankwidth,0.]],
     59            'n4':[[floodplain_width/2.+ chan_width/2.-bankwidth, floodplain_length],
     60                  [floodplain_width/2. + chan_width/2-bankwidth,0.]]
     61            }
     62
     63regionPtAreas=[[0.01, 0.01, 0.5*l0*l0],
     64               [floodplain_width/2.-chan_width/2.+0.01, 0.01, 0.5*l0*l0*0.25],
     65               [floodplain_width/2.-chan_width/2.+bankwidth+0.01, 0.01, 0.5*l0*l0],
     66               [floodplain_width/2.+chan_width/2.-bankwidth+0.01, 0.01, 0.5*l0*l0*0.25],
     67               [floodplain_width/2.+chan_width/2.+0.01, 0.01, 0.5*l0*l0]]
    6768
    6869# Define domain with appropriate boundary conditions
    69 domain = create_domain_from_regions( boundary_polygon,
    70                                    boundary_tags={'left': [0],
    71                                                   'top1': [1],
    72                                                   'chan_out': [2],
    73                                                   'top2': [3],
    74                                                   'right': [4],
    75                                                   'bottom1': [5],
    76                                                   'chan_in': [6],
    77                                                   'bottom2': [7] },
    78                                    maximum_triangle_area = 0.5*l0*l0,
    79                                    minimum_triangle_angle = 28.0,
    80                                    mesh_filename = 'channel_floodplain1.msh',
    81                                    interior_regions = [ ],
    82                                    #interior_regions = [\
    83                                    #    [channel_polygon, 0.5*l0*l0] ],
    84                                    use_cache=False,
    85                                    verbose=True)
     70anuga.create_mesh_from_regions(boundary_polygon,
     71                         boundary_tags={'left': [0],
     72                                        'top1': [1],
     73                                        'chan_out': [2],
     74                                        'top2': [3],
     75                                        'right': [4],
     76                                        'bottom1': [5],
     77                                        'chan_in': [6],
     78                                        'bottom2': [7] },
     79                         maximum_triangle_area = 0.5*l0*l0,
     80                         minimum_triangle_angle = 28.0,
     81                         filename = 'channel_floodplain1.msh',
     82                         breaklines=breakLines.values(),
     83                         regionPtArea=regionPtAreas,
     84                         use_cache=False,
     85                         verbose=True)
    8686
     87domain=anuga.create_domain_from_file('channel_floodplain1.msh')
    8788
    8889domain.set_name('channel_floodplain1') # Output name
    89 #domain.set_store_vertices_uniquely(True)
     90domain.set_store_vertices_uniquely(True)
    9091domain.set_flow_algorithm('DE1')
    91 #domain.use_edge_limiter=False
    92 #domain.extrapolate_velocity_second_order=False
    9392#------------------------------------------------------------------------------
    9493#
    … …  
    10099def topography(x, y):
    101100    # Longitudinally sloping floodplain with channel in centre
    102         #return -y*floodplain_slope -chan_bankfull_depth*\
    103     #        (x>(floodplain_width/2. - chan_width/2.) )*\
    104     #        (x<(floodplain_width/2. + chan_width/2.) )
    105 
    106101    elev1= -y*floodplain_slope - chan_bankfull_depth*\
    107102            (x>(floodplain_width/2. - chan_width/2.))*\
    … …  
    131126    return -y*floodplain_slope -chan_bankfull_depth + chan_initial_depth
    132127
    133 domain.set_quantity('elevation', topography) # Use function for elevation
     128domain.set_quantity('elevation', topography,location='centroids') # Use function for elevation
    134129domain.set_quantity('friction', man_n) # Constant friction
    135130domain.set_quantity('stage', stagetopo) # Use function for stage
    136 
    137 #def steady_uh_calc():
    138 #    # Calculate a boundary condition potentially suitable for steady state flows.
    139 #    uh = (chan_initial_depth**(4./3.)*man_n**(-2.)*floodplain_slope)**0.5\
    140 #            *chan_initial_depth  # (u) * h
    141 #    return uh
    142 #
    143 #momin = steady_uh_calc()
    144 #velin = momin/chan_initial_depth
    145 #
    146 #print 'momentum inflow is = ' + str(momin)
    147131
    148132# Define inlet operator
    … …  
    152136              [floodplain_width/2. + chan_width/2., flow_in_yval] \
    153137              ]
    154     Qin = 0.5*(floodplain_slope*(chan_width*chan_initial_depth)**2.*man_n**(-2.)\
    155             *chan_initial_depth**(4./3.) )**0.5
     138    Qin=4.5
    156139   
    157140    Inlet_operator(domain, line1, Qin)
    158141   
    159     print 'Discharge in = ', Qin #,'Velocity at inlet should be ', Qin/(chan_width*chan_initial_depth), \
    160             #'for rectangular channels, and approximately that for trapezoidal ones'
     142    print 'Discharge in = ', Qin
    161143#------------------------------------------------------------------------------
    162144#
    … …  
    165147#------------------------------------------------------------------------------
    166148
    167 # We wish to define one inflow boundary at the upstream head of the channel,
    168 # one at the outflow of the channel (e.g. transmissive)
    169 # one at the floodplain edges, excluding the downstream floodplain edge
    170 #(e.g. reflective)  and one at the downstream floodplain edge (e.g. transmissive)
    171 
    172 #Bin_sub = anuga.shallow_water.boundaries.Dirichlet_discharge_boundary(domain,\
    173 #    chan_initial_depth - chan_bankfull_depth , 0.) # An inflow for subcritical
    174 #flow
    175 #Bin_sup = anuga.shallow_water.boundaries.Dirichlet_discharge_boundary(domain, \
    176 #        chan_initial_depth - chan_bankfull_depth , steady_uh_calc()) # An inflow 
    177 #                                                                     # for
    178 #                                                                     #supercritical
    179 #                                                                     # flow
    180 #Bin_sub = anuga.Dirichlet_boundary([chan_initial_depth - chan_bankfull_depth , \
    181 #    0., 0.]) # An inflow for subcritical flow -- maybe should use transmissive
    182              # boundary instead
    183 #Bin_sup = anuga.Dirichlet_boundary([chan_initial_depth - chan_bankfull_depth ,\
    184 #        0., steady_uh_calc()]) # An inflow for supercritical flow
    185 
    186 # Trial a novel boundary, which variously imposes 2 or 3 conserved variables
    187 # depending on the local characteristic speed
    188 #Bin_trial = anuga.shallow_water.boundaries.\
    189 #        Dirichlet_inflow_boundary_sub_or_super_critical(domain, \
    190 #                stage0=chan_initial_depth - chan_bankfull_depth, wh0 = momin,\
    191 #                href = chan_initial_depth)
    192 
    193 # Trial a set_stage_transmissive_momentum inflow boundary
    194 #def inflow_stage_boundary(t):
    195 #    return chan_initial_depth - chan_bankfull_depth
    196 
    197 # transmissive_momentum_set_stage
    198 #Bin_tmss = anuga.shallow_water.boundaries.Transmissive_momentum_set_stage_boundary(domain=domain, function = inflow_stage_boundary)
    199 
    200149Br = anuga.Reflective_boundary(domain) # Solid reflective wall
    201150Bt = anuga.Transmissive_boundary(domain) # Transmissive boundary
    202 
    203 
    204151Bout_sub = anuga.Dirichlet_boundary( \
    205152        [-floodplain_length*floodplain_slope - chan_bankfull_depth + \
    … …  
    210157            + chan_initial_depth - chan_bankfull_depth
    211158
    212 # Note that the outflow boundary may be slightly incorrect for the trapezoidal channel case,
    213 # or incorrect more generally if there are numerical problems. But, in the central regions of
    214 # the channel, this shouldn't really prevent us reaching steady, uniform flow.
    215159Bout_tmss = anuga.shallow_water.boundaries.Transmissive_momentum_set_stage_boundary(domain, function = outflow_stage_boundary)
    216160
    … …  
    224168                     'chan_in': Br})
    225169
    226 
    227 # To monitor the boundary conditions, let's find centroid points representing sub and super critical flows.
    228 #inpt = numpy.argmin( ( abs(domain.centroid_coordinates[:,0] - floodplain_width/2.) + \
    229 #        abs(domain.centroid_coordinates[:,1] - 0.0)  ) )
    230 #outpt = numpy.argmin( ( abs(domain.centroid_coordinates[:, 0] - floodplain_width/2.) +\
    231 #        abs(domain.centroid_coordinates[:, 1] - floodplain_length)  ) )
    232170
    233171#------------------------------------------------------------------------------
    … …  
    257195        print '$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$'
    258196
    259     #vv = domain.get_flow_through_cross_section
    260     #print domain.quantities['ymomentum'].get_integral()
    261     #print 'Qin = ', Qin
    262     #print domain.quantities['stage'].centroid_values[inpt]
    263     #print domain.quantities['xmomentum'].centroid_values[inpt]
    264     #print domain.quantities['ymomentum'].centroid_values[inpt]
    265197
    266     # Alter inflow boundary for sub and super critical flow
    267     # FIXME: Varying the inflow boundary conditions does not work well with the modified ANUGA
    268     # (but works great with the standard ANUGA). Especially the supercritical inflow seems bad.
    269     # With the modified ANUGA, they oscillate (I wonder if the results are better with a finer grid?)
    270     #h= domain.quantities['stage'].centroid_values[inpt] - domain.quantities['elevation'].centroid_values[inpt]
    271     #uh=domain.quantities['xmomentum'].centroid_values[inpt] 
    272     #if( domain.g*h < (uh/h)**2):
    273     #    domain.set_boundary({'chan_in': Bin_sup})
    274     #else:
    275     #    domain.set_boundary({'chan_in': Bin_sub})
    276 
    277     # Alter outflow boundary for sub and super critical flow
    278     #h=domain.quantities['stage'].centroid_values[outpt] - domain.quantities['elevation'].centroid_values[outpt]
    279     #uh=domain.quantities['xmomentum'].centroid_values[outpt] 
    280     #if( domain.g*h < (uh/h)**2):
    281     #    domain.set_boundary({'chan_out': Bt})
    282     #else:
    283     #    domain.set_boundary({'chan_out': Bout_sub})
    284 
Note: See TracChangeset for help on using the changeset viewer.