source: trunk/anuga_work/development/2010-projects/anuga_1d/channel/var_depth_only2.py @ 7884

Last change on this file since 7884 was 7884, checked in by steve, 12 years ago

Moving 2010 project

File size: 2.6 KB
Line 
1import os
2import random
3from math import sqrt, pow, pi
4from channel_domain_Ab import *
5from Numeric import allclose, array, zeros, ones, Float, take, sqrt
6from config import g, epsilon
7
8
9print "Variable Width Only Test"
10
11# Define functions for initial quantities
12def initial_area(x):
13    y = zeros(len(x),Float)
14    for i in range(len(x)):
15        #if x[i]>=500:
16        #    y[i]=0
17        #else:
18       # if x[i]<900:
19        y[i]=(5-bed([x[i]]))*bed([x[i]])
20       # else:
21         #   y[i]=10-bed([x[i]])
22    return y
23
24def bed(x):
25    y = zeros(len(x),Float)
26    for i in range(len(x)):
27        if x[i]>100:
28            y[i]=2
29        else:
30            y[i]=1
31    return y
32
33def width(x):
34    return 1
35
36
37 
38import time
39
40# Set final time and yield time for simulation
41finaltime =100.0
42yieldstep = finaltime
43
44# Length of channel (m)
45L = 1000.0   
46# Define the number of cells
47number_of_cells = [10]
48
49# Define cells for finite volume and their size
50N = int(number_of_cells[0])     
51print "Evaluating domain with %d cells" %N
52cell_len = L/N # Origin = 0.0
53points = zeros(N+1,Float)
54
55# Define the centroid points
56for j in range(N+1):
57    points[j] = j*cell_len
58
59# Create domain with centroid points as defined above       
60domain = Domain(points)
61
62
63# Set initial values of quantities - default to zero
64domain.set_quantity('area', initial_area)
65domain.set_quantity('elevation',bed)
66domain.set_quantity('width',bed)
67
68# Set boundry type, order, timestepping method and limiter
69domain.set_boundary({'exterior':Reflective_boundary(domain)})
70domain.order = 2
71domain.set_timestepping_method('rk2')
72domain.set_CFL(1.0)
73domain.set_limiter("vanleer")
74#domain.h0=0.0001
75
76# Start timer
77t0 = time.time()
78
79for t in domain.evolve(yieldstep = yieldstep, finaltime = finaltime):
80    domain.write_time()
81
82N = float(N)
83HeightC = domain.quantities['height'].centroid_values
84DischargeC = domain.quantities['discharge'].centroid_values
85C = domain.centroids
86print 'That took %.2f seconds' %(time.time()-t0)
87X = domain.vertices
88HeightQ = domain.quantities['height'].vertex_values
89VelocityQ = domain.quantities['velocity'].vertex_values
90x = X.flat
91z = domain.quantities['elevation'].vertex_values.flat
92stage=HeightQ.flat+z
93b = domain.quantities['width'].vertex_values.flat
94
95
96from pylab import plot,title,xlabel,ylabel,legend,savefig,show,hold,subplot
97
98hold(False)
99   
100plot1 = subplot(211)
101
102plot(x,z,x,stage,x,b)
103 
104plot1.set_ylim([-1,11])
105xlabel('Position')
106ylabel('Stage')
107## legend(('Analytical Solution', 'Numerical Solution'),
108##            'upper right', shadow=True)
109plot2 = subplot(212)
110plot(x,VelocityQ.flat)
111plot2.set_ylim([-10,10])
112   
113xlabel('Position')
114ylabel('Velocity')
115   
116show()
Note: See TracBrowser for help on using the repository browser.