source: trunk/anuga_work/development/2010-projects/anuga_1d/sww/run_dry_dam.py @ 8073

Last change on this file since 8073 was 8073, checked in by steve, 13 years ago

Added in a unit test for inlet operator

File size: 3.1 KB
Line 
1import os
2from math import sqrt, pi
3import numpy
4import time
5#from Numeric import allclose, array, zeros, ones, Float, take, sqrt
6
7
8
9from anuga_1d.sww.sww_domain import *
10from anuga_1d.config import g, epsilon
11from anuga_1d.base.generic_mesh import uniform_mesh
12
13h1 = 10.0
14h0 = 0.1
15
16def analytical_sol(C,t):
17   
18    #t  = 0.0     # time (s)
19    # gravity (m/s^2)
20    #h1 = 10.0    # depth upstream (m)
21    #h0 = 0.0     # depth downstream (m)
22    L = 2000.0   # length of stream/domain (m)
23    n = len(C)    # number of cells
24
25    u = zeros(n,Float)
26    h = zeros(n,Float)
27    x = C-3*L/4.0
28   
29
30    for i in range(n):
31        # Calculate Analytical Solution at time t > 0
32        u3 = 2.0/3.0*(sqrt(g*h1)+x[i]/t) 
33        h3 = 4.0/(9.0*g)*(sqrt(g*h1)-x[i]/(2.0*t))*(sqrt(g*h1)-x[i]/(2.0*t))
34        u3_ = 2.0/3.0*((x[i]+L/2.0)/t-sqrt(g*h1))
35        h3_ = 1.0/(9.0*g)*((x[i]+L/2.0)/t+2*sqrt(g*h1))*((x[i]+L/2.0)/t+2*sqrt(g*h1))
36
37        if ( x[i] <= -1*L/2.0+2*(-sqrt(g*h1)*t)):
38            u[i] = 0.0
39            h[i] = h0
40        elif ( x[i] <= -1*L/2.0-(-sqrt(g*h1)*t)):
41            u[i] = u3_
42            h[i] = h3_
43
44        elif ( x[i] <= -t*sqrt(g*h1) ):
45            u[i] = 0.0 
46            h[i] = h1
47        elif ( x[i] <= 2.0*t*sqrt(g*h1) ):
48            u[i] = u3
49            h[i] = h3
50        else:
51            u[i] = 0.0 
52            h[i] = h0
53
54    return h , u*h, u
55
56
57
58def stage(x):
59    import numpy
60
61    y = numpy.where( (x>=L/4.0) & (x<=3*L/4.0), h1 , h0)
62
63#    for i in range(len(x)):
64#        if x[i]<=L/4.0:
65#            y[i] = h0
66#        elif x[i]<=3*L/4.0:
67#            y[i] = h1
68#        else:
69#            y[i] = h0
70    return y
71
72
73
74print "TEST 1D-SOLUTION III -- DRY BED"
75
76
77finaltime = 20.0
78yieldstep = 1.0
79L = 2000.0     # Length of channel (m)
80
81k = 0
82
83N = 10000
84print "Evaluating domain with %d cells" %N
85domain = Domain(*uniform_mesh(N,x_1=L))
86   
87domain.set_quantity('stage', stage)
88
89Br = Reflective_boundary(domain)
90
91domain.set_boundary({'left': Br, 'right' : Br})
92domain.set_spatial_order(2)
93domain.set_timestepping_method('rk2')
94domain.set_CFL(1.0)
95domain.set_limiter("minmod_kurganov")
96#domain.h0=0.0001
97
98t0 = time.time()
99
100for t in domain.evolve(yieldstep = yieldstep, finaltime = finaltime):
101    domain.write_time()
102
103print 'That took %.2f seconds' %(time.time()-t0)
104
105
106N = float(N)
107HeightC    = domain.quantities['height'].centroid_values
108StageC     = domain.quantities['stage'].centroid_values
109BedC       = domain.quantities['elevation'].centroid_values
110C          = domain.centroids
111
112HeightV    = domain.quantities['height'].vertex_values
113StageV     = domain.quantities['stage'].vertex_values
114BedV       = domain.quantities['elevation'].vertex_values
115VelocityV  = domain.quantities['velocity'].vertex_values
116X          = domain.vertices
117
118
119import pylab
120
121pylab.hold(False)
122
123plot1 = pylab.subplot(211)
124
125pylab.plot(X.flat,BedV.flat,X.flat,StageV.flat)
126
127plot1.set_ylim([-1,11])
128
129pylab.xlabel('Position')
130pylab.ylabel('Stage')
131pylab.legend(('Analytical Solution', 'Numerical Solution'),
132           'upper right', shadow=True)
133
134
135plot2 = pylab.subplot(212)
136
137pylab.plot(X.flat,VelocityV.flat)
138plot2.set_ylim([-15,15])
139
140pylab.xlabel('Position')
141pylab.ylabel('Velocity')
142
143pylab.show()
144
Note: See TracBrowser for help on using the repository browser.