Changeset 9223


Ignore:
Timestamp:
Jun 28, 2014, 10:30:42 AM (10 years ago)
Author:
steve
Message:

Added multiple cross sections

Location:
trunk/anuga_core
Files:
12 edited

Legend:

Unmodified
Added
Removed

  • trunk/anuga_core/source/anuga/shallow_water/sww_interrogate.py

    r8996 r9223  
    132132
    133133    return time, Q
     134
     135def get_flow_through_multiple_cross_sections(filename, polylines, verbose=False):
     136    """Obtain flow (m^3/s) perpendicular to specified cross sections.
     137
     138    filename  path to SWW file to read
     139    polylines  representation of desired cross-sections - each cross section may contain
     140              multiple sections allowing for complex shapes. Assume
     141              absolute UTM coordinates.
     142              polyline format [[x0, y0], [x1, y1], ...]
     143              polylines format [polyline_0, polyline_1, ...,polyline_n-1]
     144    verbose   True if this function is to be verbose
     145
     146    Return (time, Q)
     147    where time is a list of all stored times in SWW file
     148      and Q is a list of n hydrographs of total flow across given polylines for all
     149            stored times.
     150
     151    The normal flow is computed for each triangle intersected by the polyline
     152    and added up.  Multiple segments at different angles are specified the
     153    normal flows may partially cancel each other.
     154
     155    The typical usage of this function would be to get flow through a channel,
     156    and the polyline would then be a cross section perpendicular to the flow.
     157    """
     158
     159    quantity_names =['elevation',
     160                     'stage',
     161                     'xmomentum',
     162                     'ymomentum']
     163
     164    # Get values for quantities at each midpoint of poly line from sww file
     165    X = get_interpolated_quantities_at_multiple_polyline_midpoints(filename,
     166                                                          quantity_names=\
     167                                                              quantity_names,
     168                                                          polylines=polylines,
     169                                                          verbose=verbose)
     170    mult_segments, interpolation_function = X
     171
     172    # Get vectors for time and interpolation_points
     173    time = interpolation_function.time
     174    interpolation_points = interpolation_function.interpolation_points
     175
     176    if verbose: log.critical('Computing hydrographs')
     177
     178    # Compute hydrograph
     179    mult_Q = []
     180    base_id = 0
     181    if verbose: print '',
     182    for segments in mult_segments:
     183        if verbose: print '\b.',
     184        Q = []
     185        for t in time:
     186            total_flow = 0.0
     187            i= base_id
     188            for segment in segments:
     189                elevation, stage, uh, vh = interpolation_function(t, point_id=i)
     190                normal = segment.normal
     191
     192                # Inner product of momentum vector with segment normal [m^2/s]
     193                normal_momentum = uh*normal[0] + vh*normal[1]
     194
     195                # Flow across this segment [m^3/s]
     196                segment_flow = normal_momentum * segment.length
     197
     198                # Accumulate
     199                total_flow += segment_flow
     200               
     201                i=i+1
     202
     203            # Store flow at this timestep
     204            Q.append(total_flow)
     205           
     206        base_id = base_id + len(segments)
     207        mult_Q.append(Q)
     208
     209    if verbose: print
     210   
     211    return time, mult_Q
     212
     213def get_interpolated_quantities_at_multiple_polyline_midpoints(filename,
     214                                                      quantity_names=None,
     215                                                      polylines=None,
     216                                                      verbose=False):
     217    """Get values for quantities interpolated to multiple polyline midpoints from SWW.
     218
     219    filename        path to file to read
     220    quantity_names  quantity names to get
     221    polylines       representation of desired cross-sections
     222                    may contain multiple cross sections with multiple
     223                    sections allowing complex shapes
     224                    assume UTM coordinates
     225    verbose         True if this function is to be verbose
     226
     227    Returns (segments, i_func)
     228    where a list of segments where segments are a list of Triangle_intersection instances
     229      and i_func is an instance of Interpolation_function.
     230
     231    Note: For 'polylines' assume absolute UTM coordinates.
     232
     233    This function is used by get_flow_through_multiple_cross_sections and
     234    get_energy_through_multipe_cross_sections.
     235    """
     236
     237    from anuga.fit_interpolate.interpolate import Interpolation_function
     238
     239    # Get mesh and quantities from sww file
     240    if verbose: print 'Reading mesh and quantities from sww file'
     241    X = get_mesh_and_quantities_from_file(filename,
     242                                          quantities=quantity_names,
     243                                          verbose=verbose)
     244    mesh, quantities, time = X
     245
     246    # Find all intersections and associated triangles.
     247    mult_segments = []
     248    if verbose: print 'Intersecting segments with triangles'
     249    for polyline in polylines:
     250        mult_segments.append(mesh.get_intersecting_segments(polyline, verbose=verbose))
     251
     252    # Get midpoints
     253    interpolation_points = []
     254    for segments in mult_segments:
     255        mid_points = segment_midpoints(segments)
     256        #print mid_points
     257        interpolation_points = interpolation_points + mid_points
     258     
     259
     260    if verbose:
     261        print 'len interpolating points ', len(interpolation_points)
     262
     263    # Interpolate
     264    if verbose:
     265        log.critical('Interpolating - total number of interpolation points = %d'
     266                     % len(interpolation_points))
     267
     268    I = Interpolation_function(time,
     269                               quantities,
     270                               quantity_names=quantity_names,
     271                               vertex_coordinates=mesh.nodes,
     272                               triangles=mesh.triangles,
     273                               interpolation_points=interpolation_points,
     274                               verbose=verbose)
     275
     276    #if verbose:
     277    #    print mult_segments
     278       
     279    return mult_segments, I
     280
    134281
    135282

  • trunk/anuga_core/source/anuga/utilities/run_anuga_script.py

    r9169 r9223  
    3131            else:
    3232                cmd = 'mpirun -np %s python %s -alg %s' % (str(np), script, str(alg))
    33             print 50*'='
    34             print 'Run '+cmd
    35             print 50*'='
     33               
     34            if verbose:
     35                print 50*'='
     36                print 'Run '+cmd
     37                print 50*'='
    3638
    3739
    … …  
    4446            else:
    4547                cmd = 'python %s -alg %s' % (script, str(alg))
    46             print 50*'='
    47             print 'Run '+cmd
    48             print 50*'='
     48           
     49            if verbose:
     50                print 50*'='
     51                print 'Run '+cmd
     52                print 50*'='
     53           
    4954            os.system(cmd)
    5055            #subprocess.call([cmd], shell=True)

  • trunk/anuga_core/source/anuga/validation_utilities/produce_report.py

    r9178 r9223  
    3131
    3232    run_anuga_script('plot_results.py', args=args)
     33   
    3334    typeset_report(verbose=verbose)
    3435

  • trunk/anuga_core/source/anuga/validation_utilities/typeset_report.py

    r9164 r9223  
    1313   
    1414    import os
     15
     16    if verbose:
     17        print 50*'='
     18        print 'Running typeset_report'
     19        print 50*'='
    1520
    1621    os.system('pdflatex -shell-escape  -interaction=batchmode %s.tex' % report_name)

  • trunk/anuga_core/validation_tests/analytical_exact/runup_on_beach/numerical_runup.py

    r9155 r9223  
    99from math import sin, pi, exp
    1010from anuga import Domain
     11from anuga import myid, finalize, distribute
    1112
    12 #---------
    13 #Setup computational domain
    14 #---------
    15 points, vertices, boundary = anuga.rectangular_cross(100,3, len1=1.0,len2=0.03)
    16 domain=Domain(points,vertices,boundary)    # Create Domain
    17 domain.set_name('runup')                   # Output to file runup.sww
    18 domain.set_datadir('.')                    # Use current folder
    19 domain.set_quantities_to_be_stored({'stage': 2, 'xmomentum': 2, 'ymomentum': 2, 'elevation': 1})
    2013
    21 #------------------------------------------------------------------------------
    22 # Setup Algorithm, either using command line arguments
    23 # or override manually yourself
    24 #------------------------------------------------------------------------------
    25 from anuga.utilities.argparsing import parse_standard_args
    26 alg, cfl = parse_standard_args()
    27 domain.set_flow_algorithm(alg)
    28 #domain.set_CFL(cfl)
    29 #domain.set_minimum_allowed_height(0.0001)
    30 #domain.set_minimum_storable_height(0.0001)
     14args = anuga.get_args()
     15alg = args.alg
     16verbose = args.verbose
    3117
    32 #------------------
    33 # Define topography
    34 #------------------
    35 def topography(x,y):
    36         return -x/2 #Linear bed slope
     18if myid == 0:
     19    #---------
     20    #Setup computational domain
     21    #---------
     22    points, vertices, boundary = anuga.rectangular_cross(100,3, len1=1.0,len2=0.03)
     23    domain=Domain(points,vertices,boundary)    # Create Domain
     24    domain.set_name('runup')                   # Output to file runup.sww
     25    domain.set_datadir('.')                    # Use current folder
     26    domain.set_quantities_to_be_stored({'stage': 2, 'xmomentum': 2, 'ymomentum': 2, 'elevation': 1})
     27    domain.set_flow_algorithm(alg)
     28       
     29    #------------------
     30    # Define topography
     31    #------------------
     32    def topography(x,y):
     33            return -x/2 #Linear bed slope
     34   
     35    def stagefun(x,y):
     36        return -0.45 #Stage
     37   
     38    domain.set_quantity('elevation',topography)     # Use function for elevation
     39    domain.get_quantity('elevation').smooth_vertex_values() # Steve's fix -- without this, substantial artificial velcities are generated everywhere in the domain. With this fix, there are artificial velocities near the coast, but not elsewhere.
     40    domain.set_quantity('friction',0.0)             # Constant friction
     41    domain.set_quantity('stage', stagefun)          # Constant negative initial stage
     42else:
     43    domain = None
     44   
     45#--------------------------
     46# create Parallel Domain   
     47#--------------------------
     48domain = distribute(domain)
    3749
    38 def stagefun(x,y):
    39     return -0.45 #Stage
    40 
    41 domain.set_quantity('elevation',topography)     # Use function for elevation
    42 domain.get_quantity('elevation').smooth_vertex_values() # Steve's fix -- without this, substantial artificial velcities are generated everywhere in the domain. With this fix, there are artificial velocities near the coast, but not elsewhere.
    43 domain.set_quantity('friction',0.0)             # Constant friction
    44 domain.set_quantity('stage', stagefun)          # Constant negative initial stage
    45 
    46 #--------------------------
    4750# Setup boundary conditions
    4851#--------------------------
    … …  
    5053Bt=anuga.Transmissive_boundary(domain)          # Continue all values of boundary -- not used in this example
    5154Bd=anuga.Dirichlet_boundary([-0.4, 0., 0.])     # Constant boundary values
    52 #Bw=anuga.Time_boundary(domain=domain,
    53 #       function=lambda t: [(0.0*sin(t*2*pi)-0.1)*exp(-t)-0.1,  0.0,  0.0]) # Time varying boundary -- get rid of the 0.0 to do a runup.
    54 #        function=lambda t: [(sin(t*2*pi)-0.1)*exp(-t)-0.1,  0.0,  0.0]) # Time varying boundary
    5555
    5656#----------------------------------------------
    … …  
    7171# Produce a documentation of parameters
    7272#------------------------------------------------------------------------------
    73 parameter_file=open('parameters.tex', 'w')
    74 parameter_file.write('\\begin{verbatim}\n')
    75 from pprint import pprint
    76 pprint(domain.get_algorithm_parameters(),parameter_file,indent=4)
    77 parameter_file.write('\\end{verbatim}\n')
    78 parameter_file.close()
     73if myid == 0:
     74    parameter_file=open('parameters.tex', 'w')
     75    parameter_file.write('\\begin{verbatim}\n')
     76    from pprint import pprint
     77    pprint(domain.get_algorithm_parameters(),parameter_file,indent=4)
     78    parameter_file.write('\\end{verbatim}\n')
     79    parameter_file.close()
    7980
    8081for t in domain.evolve(yieldstep=0.2,finaltime=30.0):
    81     print domain.timestepping_statistics()
    82 print 'Finished'
     82    if myid == 0 and verbose: print domain.timestepping_statistics()
     83
     84domain.sww_merge(delete_old=True)
     85
     86finalize()
     87

  • trunk/anuga_core/validation_tests/analytical_exact/runup_on_beach/produce_results.py

    r9117 r9223  
    1 #--------------------------------
    2 # import modules
    3 #--------------------------------
    4 from anuga.validation_utilities.fabricate import *
    5 from anuga.validation_utilities import run_validation_script
    6 from anuga.validation_utilities import typeset_report
     1import anuga
     2from anuga.validation_utilities import produce_report
    73
     4args = anuga.get_args()
    85
    9 # Setup the python scripts which produce the output for this
    10 # validation test
    11 def build():
    12     run_validation_script('numerical_runup.py')
    13     run_validation_script('plot_results.py')
    14     typeset_report()
     6produce_report('numerical_runup.py', args=args)
    157
    16 def clean():
    17     autoclean()
    18 
    19 main()

  • trunk/anuga_core/validation_tests/analytical_exact/runup_on_sinusoid_beach/numerical_runup_sinusoid.py

    r9155 r9223  
    99from math import sin, pi, exp
    1010from anuga import Domain
     11from anuga import myid, finalize, distribute
    1112
    12 #---------
    13 #Setup computational domain
    14 #---------
    15 domain = anuga.rectangular_cross_domain(40,40, len1=1., len2=1.)
     13args = anuga.get_args()
     14alg = args.alg
     15verbose = args.verbose
    1616
     17scale_me=1.0
    1718
    18 #------------------------------------------------------------------------------
    19 # Setup Algorithm, either using command line arguments
    20 # or override manually yourself
    21 #------------------------------------------------------------------------------
    22 from anuga.utilities.argparsing import parse_standard_args
    23 alg, cfl = parse_standard_args()
    24 domain.set_flow_algorithm(alg)
    25 #domain.set_CFL(cfl)
    26 domain.set_name('runup_sinusoid')                         # Output to file runup.sww
    27 domain.set_datadir('.')                          # Use current folder
    28 #domain.set_quantities_to_be_stored({'stage': 2, 'xmomentum': 2, 'ymomentum': 2, 'elevation': 1})
    29 #domain.set_minimum_allowed_height(0.0001)
    30 #domain.set_minimum_storable_height(0.0001)
     19#-------------------------
     20# create Sequential domain
     21#-------------------------
     22if myid == 0:
     23        #---------
     24        #Setup computational domain
     25        #---------
     26        domain = anuga.rectangular_cross_domain(40,40, len1=1., len2=1.)
     27       
    3128
     29        domain.set_flow_algorithm(alg)
     30        domain.set_name('runup_sinusoid')                         # Output to file runup.sww
     31        domain.set_datadir('.')                          # Use current folder
     32       
     33        #------------------
     34        # Define topography
     35        #------------------
     36        def topography(x,y):
     37                return (-x/2.0 +0.05*numpy.sin((x+y)*50.0))*scale_me
     38       
     39        def stagefun(x,y):
     40                stge=-0.2*scale_me
     41                return stge
     42       
     43        domain.set_quantity('elevation',topography)     # Use function for elevation
     44        domain.get_quantity('elevation').smooth_vertex_values()
     45        domain.set_quantity('friction',0.0)             # Constant friction
     46        domain.set_quantity('stage', stagefun)          # Constant negative initial stage
     47        domain.get_quantity('stage').smooth_vertex_values()
    3248
    33 #------------------
    34 # Define topography
    35 #------------------
    36 scale_me=1.0
    37 def topography(x,y):
    38         return (-x/2.0 +0.05*numpy.sin((x+y)*50.0))*scale_me
    39 
    40 def stagefun(x,y):
    41     stge=-0.2*scale_me
    42     return stge
    43 
    44 domain.set_quantity('elevation',topography)     # Use function for elevation
    45 domain.get_quantity('elevation').smooth_vertex_values()
    46 domain.set_quantity('friction',0.0)             # Constant friction
    47 domain.set_quantity('stage', stagefun)          # Constant negative initial stage
    48 domain.get_quantity('stage').smooth_vertex_values()
    49 
     49else:
     50        domain = None
     51       
     52#------------------------
     53# Create parallel domain
     54#------------------------
     55domain = distribute(domain)
     56       
    5057
    5158#--------------------------
    … …  
    6572# Produce a documentation of parameters
    6673#------------------------------------------------------------------------------
    67 parameter_file=open('parameters.tex', 'w')
    68 parameter_file.write('\\begin{verbatim}\n')
    69 from pprint import pprint
    70 pprint(domain.get_algorithm_parameters(),parameter_file,indent=4)
    71 parameter_file.write('\\end{verbatim}\n')
    72 parameter_file.close()
     74if myid == 0:
     75        parameter_file=open('parameters.tex', 'w')
     76        parameter_file.write('\\begin{verbatim}\n')
     77        from pprint import pprint
     78        pprint(domain.get_algorithm_parameters(),parameter_file,indent=4)
     79        parameter_file.write('\\end{verbatim}\n')
     80        parameter_file.close()
    7381
    7482#------------------------------
    … …  
    7684#------------------------------
    7785for t in domain.evolve(yieldstep=1.0,finaltime=40.0):
    78     print domain.timestepping_statistics()
    79     xx = domain.quantities['xmomentum'].centroid_values
    80     yy = domain.quantities['ymomentum'].centroid_values
    81     dd = domain.quantities['stage'].centroid_values - domain.quantities['elevation'].centroid_values
    82     dd = (dd)*(dd>1.0e-03)+1.0e-06
    83     vv = ( (xx/dd)**2 + (yy/dd)**2 )**0.5
    84     vv = vv*(dd>1.0e-03)
    85     print 'Peak velocity is: ', vv.max(), vv.argmax()
    86 print 'Finished'
     86       
     87        if myid == 0 and verbose: print domain.timestepping_statistics()
     88#       xx = domain.quantities['xmomentum'].centroid_values
     89#       yy = domain.quantities['ymomentum'].centroid_values
     90#       dd = domain.quantities['stage'].centroid_values - domain.quantities['elevation'].centroid_values
     91#       dd = (dd)*(dd>1.0e-03)+1.0e-06
     92#       vv = ( (xx/dd)**2 + (yy/dd)**2 )**0.5
     93#       vv = vv*(dd>1.0e-03)
     94#       print 'Peak velocity is: ', vv.max(), vv.argmax()
     95
     96domain.sww_merge(delete_old=True)
     97
     98finalize()

  • trunk/anuga_core/validation_tests/analytical_exact/runup_on_sinusoid_beach/produce_results.py

    r9117 r9223  
    1 #--------------------------------
    2 # import modules
    3 #--------------------------------
    4 from anuga.validation_utilities.fabricate import *
    5 from anuga.validation_utilities import run_validation_script
    6 from anuga.validation_utilities import typeset_report
     1import anuga
     2from anuga.validation_utilities import produce_report
    73
    8 # Setup the python scripts which produce the output for this
    9 # validation test
    10 def build():
    11     run_validation_script('numerical_runup_sinusoid.py')
    12     run_validation_script('plot_results.py')
    13     typeset_report()
     4args = anuga.get_args()
    145
    15 def clean():
    16     autoclean()
    17 
    18 main()
     6produce_report('numerical_runup_sinusoid.py', args=args)
    197
    208
     9
     10
     11

  • trunk/anuga_core/validation_tests/analytical_exact/subcritical_over_bump/numerical_subcritical.py

    r9155 r9223  
    1010import anuga
    1111from anuga import Domain as Domain
     12from anuga import myid, finalize, distribute
    1213from math import cos
    1314from numpy import zeros, ones, float
    … …  
    2627output_file = 'subcritical'
    2728
    28 #anuga.copy_code_files(output_dir,__file__)
    29 #start_screen_catcher(output_dir+'_')
     29args = anuga.get_args()
     30alg = args.alg
     31verbose = args.verbose
    3032
    3133
    … …  
    3840W = 3*dx
    3941
    40 # structured mesh
    41 points, vertices, boundary = anuga.rectangular_cross(int(L/dx), int(W/dy), L, W, (0.0, 0.0))
    42 
    43 #domain = anuga.Domain(points, vertices, boundary)
    44 domain = Domain(points, vertices, boundary)
    45 
    46 domain.set_name(output_file)               
    47 domain.set_datadir(output_dir)
    48 
    49 #------------------------------------------------------------------------------
    50 # Setup Algorithm, either using command line arguments
    51 # or override manually yourself
    52 #------------------------------------------------------------------------------
    53 from anuga.utilities.argparsing import parse_standard_args
    54 alg, cfl = parse_standard_args()
    55 domain.set_flow_algorithm(alg)
    56 #domain.set_CFL(cfl)
    57 
    58 #------------------------------------------------------------------------------
    59 # Setup initial conditions
    60 #------------------------------------------------------------------------------
    61 def elevation(x,y):
    62     z_b = zeros(len(x))
    63     for i in range(len(x)):
    64         if (8.0 <= x[i] <= 12.0):
    65             z_b[i] = 0.2 - 0.05*(x[i]-10.0)**2.0
    66         else:
    67             z_b[i] = 0.0
    68     return z_b
    69 domain.set_quantity('elevation',elevation)
    70 domain.set_quantity('friction', 0.0)
    71 
    72 
    73 def height(x,y):
    74     return 2.0*ones(len(x))
    75 domain.set_quantity('stage', height)
     42if myid == 0:
     43    # structured mesh
     44    points, vertices, boundary = anuga.rectangular_cross(int(L/dx), int(W/dy), L, W, (0.0, 0.0))
     45   
     46    #domain = anuga.Domain(points, vertices, boundary)
     47    domain = Domain(points, vertices, boundary)
     48   
     49    domain.set_name(output_file)               
     50    domain.set_datadir(output_dir)
     51    domain.set_flow_algorithm(alg)
     52   
     53    #------------------------------------------------------------------------------
     54    # Setup initial conditions
     55    #------------------------------------------------------------------------------
     56    def elevation(x,y):
     57        z_b = zeros(len(x))
     58        for i in range(len(x)):
     59            if (8.0 <= x[i] <= 12.0):
     60                z_b[i] = 0.2 - 0.05*(x[i]-10.0)**2.0
     61            else:
     62                z_b[i] = 0.0
     63        return z_b
     64    domain.set_quantity('elevation',elevation)
     65    domain.set_quantity('friction', 0.0)
     66   
     67   
     68    def height(x,y):
     69        return 2.0*ones(len(x))
     70    domain.set_quantity('stage', height)
     71else:
     72    domain = None
     73   
     74#---------------------------
     75# Create Parallel Domain
     76#---------------------------   
     77domain = distribute(domain)
    7678
    7779#-----------------------------------------------------------------------------
    … …  
    8789
    8890
    89 #===============================================================================
    90 ##from anuga.visualiser import RealtimeVisualiser
    91 ##vis = RealtimeVisualiser(domain)
    92 ##vis.render_quantity_height("stage", zScale =h0*500, dynamic=True)
    93 ##vis.colour_height_quantity('stage', (0.0, 0.5, 1.0))
    94 ##vis.start()
    95 #===============================================================================
    9691
    9792
    … …  
    9994# Produce a documentation of parameters
    10095#------------------------------------------------------------------------------
    101 parameter_file=open('parameters.tex', 'w')
    102 parameter_file.write('\\begin{verbatim}\n')
    103 from pprint import pprint
    104 pprint(domain.get_algorithm_parameters(),parameter_file,indent=4)
    105 parameter_file.write('\\end{verbatim}\n')
    106 parameter_file.close()
     96if myid == 0:
     97    parameter_file=open('parameters.tex', 'w')
     98    parameter_file.write('\\begin{verbatim}\n')
     99    from pprint import pprint
     100    pprint(domain.get_algorithm_parameters(),parameter_file,indent=4)
     101    parameter_file.write('\\end{verbatim}\n')
     102    parameter_file.close()
    107103
    108104#------------------------------------------------------------------------------
    109105# Evolve system through time
    110106#------------------------------------------------------------------------------
    111 for t in domain.evolve(yieldstep = 1.0, finaltime = 300.):
     107for t in domain.evolve(yieldstep = 1.0, finaltime = 100.):
    112108    #print domain.timestepping_statistics(track_speeds=True)
    113     print domain.timestepping_statistics()
    114     #vis.update()
     109    if myid == 0 and verbose: print domain.timestepping_statistics()
    115110
    116111
    117 #test against know data
    118    
    119 #vis.evolveFinished()
     112domain.sww_merge(delete_old=True)
    120113
     114
     115finalize()

  • trunk/anuga_core/validation_tests/analytical_exact/subcritical_over_bump/plot_results.py

    r8797 r9223  
    66from matplotlib import pyplot as pyplot
    77from analytical_subcritical import *
    8 from numpy import ones
     8from numpy import ones, arange
    99
    1010p_st = util.get_output('subcritical.sww')
    … …  
    1212
    1313v = p2_st.y[10]
    14 v2=(p2_st.y==v)
     14#v2=(p2_st.y==v)
     15v#2=(p2_st.y>-1.0)
     16v2 = arange(len(p2_st.y))
    1517
    1618h,z = analytic_sol(p2_st.x[v2])
    … …  
    1820#Plot the stages##############################################################
    1921pyplot.clf()
    20 pyplot.plot(p2_st.x[v2], p2_st.stage[300,v2], 'b.-', label='numerical stage') # 0*T/6
     22pyplot.plot(p2_st.x[v2], p2_st.stage[-1,v2], 'b.-', label='numerical stage') # 0*T/6
    2123pyplot.plot(p2_st.x[v2], h+z,'r-', label='analytical stage')
    2224pyplot.plot(p2_st.x[v2], z,'k-', label='bed elevation')
    23 pyplot.title('Stage at an instant in time')
     25pyplot.title('Stage at time = %s secs'% p2_st.time[-1])
    2426##pyplot.ylim(-5.0,5.0)
    2527pyplot.legend(loc='best')
    … …  
    3133#Plot the momentums##########################################################
    3234pyplot.clf()
    33 pyplot.plot(p2_st.x[v2], p2_st.xmom[300,v2], 'b.-', label='numerical') # 0*T/6
     35pyplot.plot(p2_st.x[v2], p2_st.xmom[-1,v2], 'b.-', label='numerical') # 0*T/6
    3436pyplot.plot(p2_st.x[v2], 4.42*ones(len(p2_st.x[v2])),'r-', label='analytical')
    35 pyplot.title('Xmomentum at an instant in time')
     37pyplot.title('Xmomentum at time = %s secs'% p2_st.time[-1])
    3638pyplot.legend(loc='best')
    3739##pyplot.ylim([1.52,1.54])
    … …  
    4446#Plot the velocities#########################################################
    4547pyplot.clf()
    46 pyplot.plot(p2_st.x[v2], p2_st.xvel[300,v2], 'b.-', label='numerical') # 0*T/6
     48pyplot.plot(p2_st.x[v2], p2_st.xvel[-1,v2], 'b.-', label='numerical') # 0*T/6
    4749pyplot.plot(p2_st.x[v2], 4.42/h,'r-', label='analytical')
    48 pyplot.title('Xvelocity at an instant in time')
     50pyplot.title('Xvelocity at time = %s secs'% p2_st.time[-1])
    4951pyplot.legend(loc='best')
    5052pyplot.xlabel('Xposition')

  • trunk/anuga_core/validation_tests/analytical_exact/subcritical_over_bump/produce_results.py

    r9117 r9223  
    1 #--------------------------------
    2 # import modules
    3 #--------------------------------
    4 from anuga.validation_utilities.fabricate import *
    5 from anuga.validation_utilities import run_validation_script
    6 from anuga.validation_utilities import typeset_report
     1import anuga
     2from anuga.validation_utilities import produce_report
     3
     4args = anuga.get_args()
     5
     6produce_report('numerical_subcritical.py', args=args)
    77
    88
    9 # Setup the python scripts which produce the output for this
    10 # validation test
    11 def build():
    12     run_validation_script('numerical_subcritical.py')
    13     run_validation_script('plot_results.py')
    14     typeset_report()
    15 
    16 def clean():
    17     autoclean()
    18 
    19 main()
    20 

  • trunk/anuga_core/validation_tests/analytical_exact/trapezoidal_channel/plot_results.py

    r9221 r9223  
    7777
    7878analytical_stage = min(p.elev[v1]) + dc_analytical
    79 analytic_vel = ( (1./300.)*(analytical_stage-p.elev[v1])**(4./3.)*(1./0.03)**2.)**0.5
     79analytic_vel = ( (1./300.)*numpy.maximum(analytical_stage-p.elev[v1],0.0)**(4./3.)*(1./0.03)**2.)**0.5
    8080analytic_vel = analytic_vel*(analytical_stage>p.elev[v1])
    8181
    … …  
    107107
    108108analytical_stage = min(p.elev[v1]) + dc_analytical
    109 analytic_vel = ( (1./300.)*(analytical_stage-p.elev[v1])**(4./3.)*(1./0.03)**2.)**0.5
     109analytic_vel = ( (1./300.)*numpy.maximum(analytical_stage-p.elev[v1],0.0)**(4./3.)*(1./0.03)**2.)**0.5
    110110analytic_vel = analytic_vel*(analytical_stage>p.elev[v1])
    111111
    … …  
    138138
    139139analytical_stage = min(p.elev[v1]) + dc_analytical
    140 analytic_vel = ( (1./300.)*(analytical_stage-p.elev[v1])**(4./3.)*(1./0.03)**2.)**0.5
     140analytic_vel = ( (1./300.)*numpy.maximum(analytical_stage-p.elev[v1],0.0)**(4./3.)*(1./0.03)**2.)**0.5
    141141analytic_vel = analytic_vel*(analytical_stage>p.elev[v1])
    142142
    … …  
    167167print '#======================================================================'
    168168
    169 from anuga.shallow_water.sww_interrogate import get_flow_through_cross_section
     169from anuga.shallow_water.sww_interrogate import get_flow_through_multiple_cross_sections
    170170
    171171polyline0 = [ [floodplain_width, 10.0], [0., 10.0]]
    … …  
    173173polyline2 = [[floodplain_width, floodplain_length-1.0], [0., floodplain_length-1.0]]
    174174       
    175 time,Q0  = get_flow_through_cross_section(filename, polyline0, verbose=True)
    176 time,Q1  = get_flow_through_cross_section(filename, polyline1, verbose=True)
    177 time,Q2  = get_flow_through_cross_section(filename, polyline2, verbose=True)
     175polylines= [polyline0, polyline1, polyline2]       
     176
     177time, [Q0,Q1,Q2]  = get_flow_through_multiple_cross_sections(filename, polylines, verbose=True)
    178178
    179179pyplot.figure(figsize=(12.,8.))
Note: See TracChangeset for help on using the changeset viewer.