source: anuga_work/development/anuga_1d/analytic_dam_sudi.py @ 5535

Last change on this file since 5535 was 5535, checked in by steve, 16 years ago

Changed name of shallow_water_1d to anuga_1d

File size: 3.8 KB
Line 
1
2
3class AnalyticDam:
4
5    def __init__(self, h0 = 5.0, h1 = 10.0, L = 2000.0):
6        from math import sqrt
7       
8        self.h0 = h0 # depth upstream (m)
9        self.h1 = h1 # depth downstream (m)
10        self.= L  # length of domain
11
12        g  = 9.81    # gravity (m/s^2)
13       
14        c0 = sqrt(g*h0) #left celerity
15        c1 = sqrt(g*h1) #right celerity
16       
17        zmin=-100.0
18        zmax=101.0
19        for i in range(100):
20            z=(zmin+zmax)/2.0
21            u2=z-c0*c0/4.0/z*(1.0+sqrt(1.0+8.0*z*z/c0/c0))
22            c2=c0*sqrt(0.5*(sqrt(1.0+8.0*z*z/c0/c0)-1.0))
23            func=2.0*c1/c0-u2/c0-2.0*c2/c0
24            if (func > 0.0):
25                zmin=z
26            else:
27                zmax=z
28
29        if( abs(z) > 99.0):
30            print 'no convergence'
31
32        self.u2 = u2
33        self.c0 = c0
34        self.c1 = c1
35        self.c2 = c2
36        self.g = g
37        self.z = z
38
39
40
41    def __call__(self, C,t):
42       
43        from Numeric import zeros,Float
44        from math import sqrt
45   
46        #t  = 0.0     # time (s)
47        h0 = self.h0   
48        h1 = self.h1   
49        L = self.L
50        n = len(C)    # number of cells
51
52        u2 = self.u2
53        c0 = self.c0
54        c1 = self.c1
55        c2 = self.c2
56       
57        g = self.g
58        z = self.z
59
60        u = zeros(n,Float)
61        h = zeros(n,Float)
62        uh = zeros(n,Float)
63        x = C-3*L/4.0
64        #x = zeros(n,Float)
65        #for i in range(n):
66        #    x[i] = C[i]-1000.0
67
68        # Upstream and downstream boundary conditions are set to the intial water
69        # depth for all time.
70
71        # Calculate Shock Speed
72        #h2 = 7.2692044
73       
74        #S2 = 2*h2/(h2-h0)*(sqrt(g*h1)-sqrt(g*h2))
75        #u2 = S2 - g*h0/(4*S2)*(1+sqrt(1+8*S2*S2/(g*h0)))
76
77        h2=h0/(1.0-u2/z)
78        x3=(u2-c2)*t
79        x2=z*t
80        x1=-c1*t
81        x1_ = -1*L/2.0-x1
82        x2_ = -1*L/2.0+2*x1
83        #x3_ = -1*L/2.0-x3
84        #t=50
85        #x = (-L/2:L/2)
86        for i in range(n):
87            # Calculate Analytical Solution at time t > 0
88            u3 = 2.0/3.0*(sqrt(g*h1)+x[i]/t) 
89            h3 = 4.0/(9.0*g)*(sqrt(g*h1)-x[i]/(2.0*t))*(sqrt(g*h1)-x[i]/(2.0*t))
90            u3_ = 2.0/3.0*((x[i]+L/2.0)/t-sqrt(g*h1))
91            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))
92            #if t == 30:
93            #    x[i] = 500
94            #    print 'x2',x2
95            #    print 'x3',x3
96            #    print 'x1',x1
97            if ( x[i] <= x2_ ):
98                #print 'here x2_=', x2_
99                u[i] = 0.0
100                h[i] = 0.0
101                uh [i] = u[i]*h[i]
102            #elif ( x[i] <= x3_ ):
103            #    print 'here x3_=', x3_
104            #    u[i] = -1*u2
105            #    h[i] = h2
106            #    uh[i] = u[i]*h[i]
107            elif ( x[i] <= x1_ ):
108                #print 'here x1_=', x1_
109                u[i] = u3_
110                h[i] = h3_
111                uh[i] = u[i]*h[i]
112            #else:
113            #    u[i] = 0.0
114            #    h[i] = h0
115            #    uh[i] = u[i]*h[i]
116
117            #elif ( x[i] <= x1/2.0 ):
118            #    u[i] = 0.0
119            #    h[i] = h1
120            #    uh[i] = u[i]*h[i]
121            elif ( x[i] <= x1 ):
122                #print 'here x1=', x1
123                u[i] = 0.0 
124                h[i] = h1
125                uh[i] = u[i]*h[i]
126            elif ( x[i] <= x3 ):
127                #print 'here x3=', x3
128                u[i] = u3
129                h[i] = h3
130                uh[i] = u[i]*h[i]
131            elif ( x[i] < x2 ):
132                #print 'here x2=', x2
133                u[i] = u2
134                h[i] = h2
135                uh[i] = u[i]*h[i]
136            else:
137                #print 'here the last section'
138                u[i] = 0.0 
139                h[i] = h0
140                uh[i] = u[i]*h[i]
141
142        return h , uh
Note: See TracBrowser for help on using the repository browser.