 Timestamp:
 Jun 15, 2010, 5:34:24 PM (13 years ago)
 File:

 1 edited
Legend:
 Unmodified
 Added
 Removed

anuga_work/development/2010projects/anuga_1d/channel/channel_flow_1_padarn.py
r7839 r7842 2 2 from math import sqrt, pow, pi 3 3 4 from Numeric import allclose, array, zeros, ones, Float, take, sqrt 4 import numpy 5 5 6 from channel_domain import * 7 from config import g, epsilon 6 from anuga_1d.channel.channel_domain import * 7 from anuga_1d.config import g, epsilon 8 from anuga_1d.generic.generic_mesh import interval_mesh 8 9 9 10 … … 21 22 22 23 def 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 30 28 31 29 def initial_discharge(x): … … 35 33 36 34 # Set final time and yield time for simulation 37 finaltime = 10.038 yieldstep = finaltime35 finaltime = 50.0 36 yieldstep = 10.0 39 37 40 38 # Length of channel (m) 41 39 L = 1000.0 42 40 # 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 41 N = 200 54 42 55 43 # Create domain with centroid points as defined above 56 domain = Domain(points) 44 domain = Domain(*interval_mesh(N)) 45 57 46 58 47 # Set initial values of quantities  default to zero … … 63 52 64 53 # Set boundry type, order, timestepping method and limiter 65 domain.set_boundary({'exterior': Dirichlet_boundary([14,20,0,1.4,20/14,9])}) 54 Bd = Dirichlet_boundary([14,20,0,1.4,20/14,9,1.4]) 55 domain.set_boundary({'left': Bd , 'right' : Bd }) 56 57 66 58 domain.order = 2 67 59 domain.set_timestepping_method('rk2') … … 76 68 domain.write_time() 77 69 78 N = float(N) 79 HeightC = domain.quantities['height'].centroid_values 70 71 exit() 72 73 HeightC = domain.quantities['height'].centroid_values 80 74 DischargeC = domain.quantities['discharge'].centroid_values 81 C = domain.centroids75 C = domain.centroids 82 76 print '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 77 X = domain.vertices 78 HeightQ = domain.quantities['height'].vertex_values 79 VelocityQ = domain.quantities['velocity'].vertex_values 80 Z = domain.quantities['elevation'].vertex_values 81 Stage = HeightQ+Z 89 82 90 83 from pylab import plot,title,xlabel,ylabel,legend,savefig,show,hold,subplot … … 94 87 plot1 = subplot(211) 95 88 96 plot( x,z,x,stage)89 plot(X.flat,Z.flat,X.flat,Stage.flat) 97 90 98 91 plot1.set_ylim([1,11]) … … 101 94 legend(('Analytical Solution', 'Numerical Solution'), 102 95 'upper right', shadow=True) 96 103 97 plot2 = subplot(212) 104 plot( x,VelocityQ.flat)98 plot(X.flat,VelocityQ.flat) 105 99 plot2.set_ylim([10,10]) 106 100
Note: See TracChangeset
for help on using the changeset viewer.