Changeset 4703


Ignore:
Timestamp:
Sep 5, 2007, 2:23:41 PM (17 years ago)
Author:
sexton
Message:

updated scripts for convergence studu (note, png files created from animatesww2d_alt.py then need to be dropped into something like Windows Movie Maker to export to a mpeg)

Location:
anuga_validation/convergence_study
Files:
2 edited

Legend:

Unmodified
Added
Removed

  • anuga_validation/convergence_study/animatesww2d_alt.py

    r4700 r4703  
    1717    ymomentum = fid.variables['ymomentum']
    1818   
    19 
    20     # find indices for where y = 0
    21 
    22     #
    23     # Set up plot
    24 
    25     spacing=10 # grid interval in metres
    26     hmax=5     # ??
    27 
    28    
    29 #figure(1);
    30 
    31 #if nargin > 1
    32 #    figure(1);set(gcf,'Position',[200 400 800 300]);
    33 #    mov=avifile(movie_filename,'FPS',5)
    34 #end
    35 
    36 #i=find(V_y==0);
    37 
    38 #[xi,yi]=meshgrid(min(V_x):spacing:max(V_x),min(V_y):spacing:max(V_y));
    39 
    40 #elev=griddata(V_x,V_y,V_elevation,xi,yi);
    41 
    42 #if nargin < 3
    43 #    range=[min(V_x) max(V_x) min(min(V_stage))-1 max(max(V_stage))+1];
    44 #end
    45 
    46 #m=V_stage(1,i);
    47 
    48 #for t=1:length(V_time)
    49    
    50 #    m=max(m,V_stage(t,i));
    51 #    plot(V_x(i),V_stage(t,i),V_x(i),V_elevation(i),V_x(i),m,'.');
    52 
    53 #    axis(range)
    54    
    55 
    56 #    title(num2str(V_time(t)));drawnow;
    57 
    58 #    if nargin>1
    59 #        F=getframe(gcf);
    60 #        mov=addframe(mov,F);   
    61 #    end
    62    
    63 #end
    64 
    65 #if nargin > 1
    66 #    mov=close(mov);
    67 #end
    68 
    69 '''
    70 if __name__ == '__main__':
    71 
    72     import sys
    73     animatesww2d_alt(sys.argv[1], 'cross_section_plot', 5)
    74 
    75 '''
    76 
    77 '''
    78 How do I make a movie with matplotlib?
    79 If you want to take an animated plot and turn it into a movie, the best approach is to save a series of image files (eg PNG) and use an external tool to convert them to a movie. There is a matplotlib tutorial on this subject here. You can use mencoder, which is part of the mplayer suite for this
    80 #fps (frames per second) controls the play speed
    81 mencoder 'mf://*.png' -mf type=png:fps=10 -ovc \
    82    lavc -lavcopts vcodec=wmv2 -oac copy -o animation.avi
    83 
    84 The swiss army knife of image tools, ImageMagick's convert, works for this as well.
    85 Here is a simple example script that saves some PNGs, makes them into a movie, and then cleans up.
    86 
    87 import os, sys
    88 from pylab import *
    89 
    90 files = []
    91 figure(figsize=(5,5))
    92 ax = subplot(111)
    93 for i in range(50):  # 50 frames
    94     cla()
    95     imshow(rand(5,5), interpolation='nearest')
    96     fname = '_tmp%03d.png'%i
    97     print 'Saving frame', fname
    98     savefig(fname)
    99     files.append(fname)
    100 
    101 print 'Making movie animation.mpg - this make take a while'
    102 os.system("mencoder 'mf://_tmp*.png' -mf type=png:fps=10 \
    103   -ovc lavc -lavcopts vcodec=wmv2 -oac copy -o animation.mpg")
    104 
    105 # cleanup
    106 for fname in files: os.remove(fname)
    107 
    108 '''
    109 
    11019def animatesww2d(swwfile):
    11120    """Read in sww file and plot cross section of model output
    … …  
    12433
    12534    # interpolation points are for y = 0
    126     m = 100 # number of increments in x
     35    # number of increments in x
     36    m = 100
    12737    max_x = 100000.
    12838    x = 0
    … …  
    13141    for i in range(m):
    13242        x += max_x/m
    133         points.append([x,0])
     43        points.append([x,0.])
    13444        x_points.append(x)
    13545           
    … …  
    15767    momenta = zeros((n,m), Float)
    15868    speed = zeros((n,m), Float)
    159     max_stages = []
     69    max_stages = zeros((n,m), Float)
    16070    max_stage = None
     71    min_stages = zeros((n,m), Float)
     72    min_stage = 100.
    16173
    16274    for k, location in enumerate(x_points):
    … …  
    18193
    18294            if w > max_stage: max_stage = w
    183             max_stages.append(max_stage)
     95            if w < min_stage: min_stage = w
     96            max_stages[i,k] = max_stage
     97            min_stages[i,k] = min_stage
    18498
    185     from pylab import plot, xlabel, ylabel, savefig, close, hold
     99    min_stage = min(min(min_stages))
     100    max_stage = max(max(max_stages))
     101    stage_axis = [0, max_x, min_stage-1, max_stage+1]
     102    from pylab import plot, xlabel, ylabel, savefig, close, hold, axis, title
    186103    hold(False)
    187     for i in range(n):
    188         plot(x_points,stage[i,:])
    189         xlabel('x')
     104    figs = []
     105    j = 1
     106    #max_vec = zeros(100, Float)
     107
     108    for i in range(0,n,20):
     109        #max_vec[i] = max(stage[i,:])
     110        plot(x_points,stage[i,:]) #x_points,max_vec,'+')
     111        xlabel('x (m)')
    190112        ylabel('stage (m)')
    191         #title('')
     113        axis(stage_axis)
    192114        name = 'time_%i' %i
    193115        savefig(name)
     116        title(name)
     117        figs.append(name)
    194118
    195119    close('all')
    196120   
    197121animatesww2d('myexample2.sww')
     122
     123# once png files have been created, they can then be dragged into
     124# Windows Movie Maker and exported to a mpeg.

  • anuga_validation/convergence_study/convergence.py

    r4694 r4703  
    44similar to a beach environment
    55"""
    6 
    76
    87#------------------------------------------------------------------------------
    … …  
    1817from anuga.shallow_water import Transmissive_boundary
    1918from anuga.shallow_water import Transmissive_Momentum_Set_Stage_boundary
    20 #from anuga.fit_interpolate.fit import fit_to_mesh_file
    2119from anuga.geospatial_data.geospatial_data import *
    2220from math import cos
    … …  
    2523# Setup computational domain
    2624#------------------------------------------------------------------------------
    27 
    28 dx = 1000.
     25dx = 10.
    2926dy = dx
    3027L = 100000.
    3128W = 3000.
    3229
    33 points, vertices, boundary = rectangular_cross(int(L/dx), int(L/dy),
    34                                                L, W, (0.0, -W/2)) # Basic mesh
    35 
     30# structured mesh
     31#points, vertices, boundary = rectangular_cross(int(L/dx), int(W/dy),
     32#                                               L, W, (0.0, -W/2)) # Basic mesh
    3633#points, vertices, boundary = rectangular_cross(666, 3, 100000, 3000, (0.0, -0.0)) # Basic mesh
    3734#points, vertices, boundary = rectangular_cross(530, 10, 5300, 100, (-5000.0, -50.0)) # Basic mesh
    3835#points, vertices, boundary = rectangular_cross(1000, 100, 20, 3) # Basic mesh
     36#domain = Domain(points, vertices, boundary)
    3937
    40 domain = Domain(points, vertices, boundary) # Create domain
    41 domain.set_name('myexample2')                    # Output to bedslope.sww
     38# unstructured mesh
     39poly_domain = [[0,-W],[0,W],[L,W],[L,-W]]
     40meshname = 'test.msh'
     41from anuga.pmesh.mesh_interface import create_mesh_from_regions
     42# Create mesh
     43create_mesh_from_regions(poly_domain,
     44                         boundary_tags={'left': [0], 'top': [1],
     45                                        'right': [2], 'bottom': [3]},
     46                         maximum_triangle_area = 25000,
     47                         filename=meshname)
     48
     49# Create domain
     50domain = Domain(meshname, use_cache=True, verbose = True)
     51domain.set_name('myexample2')               
    4252domain.set_datadir('.')                     # Use current directory for output
    43 #domain.set_quantities_to_be_stored(['stage',# Store all conserved quantities
    44 #                                    'xmomentum',
    45 #                                    'ymomentum'])   
    46 
    4753
    4854#------------------------------------------------------------------------------
    4955# Setup initial conditions
    5056#------------------------------------------------------------------------------
    51 
    52 
    5357#domain.set_quantity('elevation', topography) # Use function for elevation
    5458domain.set_quantity('elevation',-100)
    55 #domain.set_quantity('elevation',filename='tr10_mod.pts')
    56 #domain.set_quantity('elevation',filename='tr5_mod.pts')
    57 #domain.set_quantity('elevation',filename='tr4_mod.pts')
    58 #domain.set_quantity('friction', 0.1)         # Constant friction
    59 domain.set_quantity('friction', 0.00)         # Constant friction
    60 #domain.set_quantity('friction', 0.035)         # Constant friction
    61 domain.set_quantity('stage', 0.0)            # Constant negative initial stage
    62 
     59domain.set_quantity('friction', 0.00)
     60domain.set_quantity('stage', 0.0)           
    6361
    6462#-----------------------------------------------------------------------------
    6563# Setup boundary conditions
    6664#------------------------------------------------------------------------------
    67 
    6865from math import sin, pi, exp
    6966Br = Reflective_boundary(domain)      # Solid reflective wall
    7067Bt = Transmissive_boundary(domain)    # Continue all values on boundary
    7168Bd = Dirichlet_boundary([1,0.,0.]) # Constant boundary values
    72 amplitude=eval(sys.argv[1])
     69amplitude = 1
    7370Bw = Time_boundary(domain=domain,     # Time dependent boundary 
    7471## Sine wave
    … …  
    8178#                   f=lambda t: [(-8.0*sin((1./720.)*t*2*pi))*((t<720.)-0.5*(t<360.)), 0.0, 0.0])
    8279
    83 
    8480# Associate boundary tags with boundary objects
    85 #domain.set_boundary({'left': Bw, 'right': Br, 'top': Br, 'bottom': Br})
    8681domain.set_boundary({'left': Bw, 'right': Bt, 'top': Br, 'bottom': Br})
    87 
    8882
    8983#------------------------------------------------------------------------------
    … …  
    9488    domain.write_time()
    9589   
    96 
Note: See TracChangeset for help on using the changeset viewer.