Ignore:
Timestamp:
Jun 15, 2010, 5:34:24 PM (13 years ago)
Author:
steve
Message:

Adding in a few new files

File:
1 edited

Legend:

Unmodified
Added
Removed
  • anuga_work/development/2010-projects/anuga_1d/channel/channel_flow_1_padarn.py

    r7839 r7842  
    22from math import sqrt, pow, pi
    33
    4 from Numeric import allclose, array, zeros, ones, Float, take, sqrt
     4import numpy
    55
    6 from channel_domain import *
    7 from config import g, epsilon
     6from anuga_1d.channel.channel_domain import *
     7from anuga_1d.config import g, epsilon
     8from anuga_1d.generic.generic_mesh import interval_mesh
    89
    910
     
    2122
    2223def bed(x):
    23     y = zeros(len(x),Float)
    24     for i in range(len(x)):
    25         if x[i]<525 and x[i]>475:
    26             y[i]=1
    27         else:
    28             y[i]=0
    29     return y
     24    y = numpy.ones(len(x),numpy.float)
     25
     26    return numpy.where( (x<525) & (x>475),y,0.0)
     27   
    3028
    3129def initial_discharge(x):
     
    3533
    3634# Set final time and yield time for simulation
    37 finaltime = 10.0
    38 yieldstep = finaltime
     35finaltime = 50.0
     36yieldstep = 10.0
    3937
    4038# Length of channel (m)
    4139L = 1000.0   
    4240# Define the number of cells
    43 number_of_cells = [200]
    44 
    45 # Define cells for finite volume and their size
    46 N = int(number_of_cells[0])     
    47 print "Evaluating domain with %d cells" %N
    48 cell_len = L/N # Origin = 0.0
    49 points = zeros(N+1,Float)
    50 
    51 # Define the centroid points
    52 for j in range(N+1):
    53     points[j] = j*cell_len
     41N = 200
    5442
    5543# Create domain with centroid points as defined above       
    56 domain = Domain(points)
     44domain = Domain(*interval_mesh(N))
     45
    5746
    5847# Set initial values of quantities - default to zero
     
    6352
    6453# Set boundry type, order, timestepping method and limiter
    65 domain.set_boundary({'exterior': Dirichlet_boundary([14,20,0,1.4,20/14,9])})
     54Bd = Dirichlet_boundary([14,20,0,1.4,20/14,9,1.4])
     55domain.set_boundary({'left': Bd , 'right' : Bd })
     56
     57
    6658domain.order = 2
    6759domain.set_timestepping_method('rk2')
     
    7668    domain.write_time()
    7769
    78 N = float(N)
    79 HeightC = domain.quantities['height'].centroid_values
     70
     71exit()
     72
     73HeightC    = domain.quantities['height'].centroid_values
    8074DischargeC = domain.quantities['discharge'].centroid_values
    81 C = domain.centroids
     75C          = domain.centroids
    8276print 'That took %.2f seconds' %(time.time()-t0)
    83 X = domain.vertices
    84 HeightQ = domain.quantities['height'].vertex_values
    85 VelocityQ = domain.quantities['velocity'].vertex_values
    86 x = X.flat
    87 z = domain.quantities['elevation'].vertex_values.flat
    88 stage=HeightQ.flat+z
     77X          = domain.vertices
     78HeightQ    = domain.quantities['height'].vertex_values
     79VelocityQ  = domain.quantities['velocity'].vertex_values
     80Z = domain.quantities['elevation'].vertex_values
     81Stage      = HeightQ+Z
    8982
    9083from pylab import plot,title,xlabel,ylabel,legend,savefig,show,hold,subplot
     
    9487plot1 = subplot(211)
    9588
    96 plot(x,z,x,stage)
     89plot(X.flat,Z.flat,X.flat,Stage.flat)
    9790 
    9891plot1.set_ylim([-1,11])
     
    10194legend(('Analytical Solution', 'Numerical Solution'),
    10295           'upper right', shadow=True)
     96
    10397plot2 = subplot(212)
    104 plot(x,VelocityQ.flat)
     98plot(X.flat,VelocityQ.flat)
    10599plot2.set_ylim([-10,10])
    106100   
Note: See TracChangeset for help on using the changeset viewer.