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

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

Added in the well balancing of the swb_domain by combining all the pressure terms

File size: 3.2 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.0
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 = 2000
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")
96domain.set_beta(1.5)
97#domain.h0=0.0001
98
99t0 = time.time()
100
101for t in domain.evolve(yieldstep = yieldstep, finaltime = finaltime):
102    domain.write_time()
103
104print 'That took %.2f seconds' %(time.time()-t0)
105
106
107N = float(N)
108HeightC    = domain.quantities['height'].centroid_values
109StageC     = domain.quantities['stage'].centroid_values
110BedC       = domain.quantities['elevation'].centroid_values
111C          = domain.centroids
112
113HeightV    = domain.quantities['height'].vertex_values
114StageV     = domain.quantities['stage'].vertex_values
115BedV       = domain.quantities['elevation'].vertex_values
116VelocityV  = domain.quantities['velocity'].vertex_values
117X          = domain.vertices
118
119
120import pylab
121
122pylab.hold(False)
123
124plot1 = pylab.subplot(211)
125
126pylab.plot(X.flat,BedV.flat,X.flat,StageV.flat)
127
128plot1.set_ylim([-1,11])
129
130pylab.xlabel('Position')
131pylab.ylabel('Stage')
132pylab.legend(('Analytical Solution', 'Numerical Solution'),
133           'upper right', shadow=True)
134
135
136plot2 = pylab.subplot(212)
137
138pylab.plot(X.flat,VelocityV.flat)
139plot2.set_ylim([-20,20])
140
141pylab.xlabel('Position')
142pylab.ylabel('Velocity')
143
144pylab.show()
145
Note: See TracBrowser for help on using the repository browser.