Changeset 7842 for anuga_work/development/2010-projects/anuga_1d/channel/channel_flow_1_padarn.py

Timestamp:

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

Author:

steve

Message:

Adding in a few new files

File:

• anuga_work/development/2010-projects/anuga_1d/channel/channel_flow_1_padarn.py (modified) (7 diffs)

anuga_work/development/2010-projects/anuga_1d/channel/channel_flow_1_padarn.py

@@ -2 +2 @@
 from math import sqrt, pow, pi
 
-from Numeric import allclose, array, zeros, ones, Float, take, sqrt
+import numpy
 
-from channel_domain import *
-from config import g, epsilon
+from anuga_1d.channel.channel_domain import *
+from anuga_1d.config import g, epsilon
+from anuga_1d.generic.generic_mesh import interval_mesh
 
 
@@ -21 +22 @@
 
 def bed(x):
-    y = zeros(len(x),Float)
-    for i in range(len(x)):
-        if x[i]<525 and x[i]>475:
-            y[i]=1
-        else:
-            y[i]=0
-    return y
+    y = numpy.ones(len(x),numpy.float)
+
+    return numpy.where( (x<525) & (x>475),y,0.0)
+    
 
 def initial_discharge(x):
@@ -35 +33 @@
 
 # Set final time and yield time for simulation
-finaltime = 10.0
-yieldstep = finaltime
+finaltime = 50.0
+yieldstep = 10.0
 
 # Length of channel (m)
 L = 1000.0   
 # Define the number of cells
-number_of_cells = [200]
-
-# Define cells for finite volume and their size
-N = int(number_of_cells[0])     
-print "Evaluating domain with %d cells" %N
-cell_len = L/N # Origin = 0.0
-points = zeros(N+1,Float)
-
-# Define the centroid points
-for j in range(N+1):
-    points[j] = j*cell_len
+N = 200
 
 # Create domain with centroid points as defined above        
-domain = Domain(points)
+domain = Domain(*interval_mesh(N))
+
 
 # Set initial values of quantities - default to zero
@@ -63 +52 @@
 
 # Set boundry type, order, timestepping method and limiter
-domain.set_boundary({'exterior': Dirichlet_boundary([14,20,0,1.4,20/14,9])})
+Bd = Dirichlet_boundary([14,20,0,1.4,20/14,9,1.4])
+domain.set_boundary({'left': Bd , 'right' : Bd })
+
+
 domain.order = 2
 domain.set_timestepping_method('rk2')
@@ -76 +68 @@
     domain.write_time()
 
-N = float(N)
-HeightC = domain.quantities['height'].centroid_values
+
+exit()
+
+HeightC    = domain.quantities['height'].centroid_values
 DischargeC = domain.quantities['discharge'].centroid_values
-C = domain.centroids
+C          = domain.centroids
 print 'That took %.2f seconds' %(time.time()-t0)
-X = domain.vertices
-HeightQ = domain.quantities['height'].vertex_values
-VelocityQ = domain.quantities['velocity'].vertex_values
-x = X.flat
-z = domain.quantities['elevation'].vertex_values.flat
-stage=HeightQ.flat+z
+X          = domain.vertices
+HeightQ    = domain.quantities['height'].vertex_values
+VelocityQ  = domain.quantities['velocity'].vertex_values
+Z = domain.quantities['elevation'].vertex_values
+Stage      = HeightQ+Z
 
 from pylab import plot,title,xlabel,ylabel,legend,savefig,show,hold,subplot
@@ -94 +87 @@
 plot1 = subplot(211)
 
-plot(x,z,x,stage)
+plot(X.flat,Z.flat,X.flat,Stage.flat)
  
 plot1.set_ylim([-1,11])
@@ -101 +94 @@
 legend(('Analytical Solution', 'Numerical Solution'),
            'upper right', shadow=True)
+
 plot2 = subplot(212)
-plot(x,VelocityQ.flat)
+plot(X.flat,VelocityQ.flat)
 plot2.set_ylim([-10,10])