source: trunk/anuga_work/development/sudi/sw_1d/shock_detector/shm/analytic_dam_sudi.py @ 7909

Last change on this file since 7909 was 7909, checked in by mungkasi, 14 years ago

Add 1d codes.

File size: 3.9 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            #print 'func=',func
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 Float
44        from numpy import zeros       
45        from math import sqrt
46   
47        #t  = 0.0     # time (s)
48        h0 = self.h0   
49        h1 = self.h1   
50        L = self.L
51        n = len(C)    # number of cells
52
53        u2 = self.u2
54        c0 = self.c0
55        c1 = self.c1
56        c2 = self.c2
57       
58        g = self.g
59        z = self.z
60
61        u = zeros(n,Float)
62        h = zeros(n,Float)
63        uh = zeros(n,Float)
64        x = C-3*L/4.0
65        #x = zeros(n,Float)
66        #for i in range(n):
67        #    x[i] = C[i]-1000.0
68
69        # Upstream and downstream boundary conditions are set to the intial water
70        # depth for all time.
71
72        # Calculate Shock Speed
73        #h2 = 7.2692044
74       
75        #S2 = 2*h2/(h2-h0)*(sqrt(g*h1)-sqrt(g*h2))
76        #u2 = S2 - g*h0/(4*S2)*(1+sqrt(1+8*S2*S2/(g*h0)))
77
78        h2=h0/(1.0-u2/z)
79        x3=(u2-c2)*t
80        x2=z*t
81        x1=-c1*t
82        x1_ = -1*L/2.0-x1
83        x2_ = -1*L/2.0+2*x1
84        #x3_ = -1*L/2.0-x3
85        #t=50
86        #x = (-L/2:L/2)
87        for i in range(n):
88            # Calculate Analytical Solution at time t > 0
89            u3 = 2.0/3.0*(sqrt(g*h1)+x[i]/t) 
90            h3 = 4.0/(9.0*g)*(sqrt(g*h1)-x[i]/(2.0*t))*(sqrt(g*h1)-x[i]/(2.0*t))
91            u3_ = 2.0/3.0*((x[i]+L/2.0)/t-sqrt(g*h1))
92            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))
93            #if t == 30:
94            #    x[i] = 500
95            #    print 'x2',x2
96            #    print 'x3',x3
97            #    print 'x1',x1
98            if ( x[i] <= x2_ ):
99                #print 'here x2_=', x2_
100                u[i] = 0.0
101                h[i] = 0.0
102                uh [i] = u[i]*h[i]
103            #elif ( x[i] <= x3_ ):
104            #    print 'here x3_=', x3_
105            #    u[i] = -1*u2
106            #    h[i] = h2
107            #    uh[i] = u[i]*h[i]
108            elif ( x[i] <= x1_ ):
109                #print 'here x1_=', x1_
110                u[i] = u3_
111                h[i] = h3_
112                uh[i] = u[i]*h[i]
113            #else:
114            #    u[i] = 0.0
115            #    h[i] = h0
116            #    uh[i] = u[i]*h[i]
117
118            #elif ( x[i] <= x1/2.0 ):
119            #    u[i] = 0.0
120            #    h[i] = h1
121            #    uh[i] = u[i]*h[i]
122            elif ( x[i] <= x1 ):
123                #print 'here x1=', x1
124                u[i] = 0.0 
125                h[i] = h1
126                uh[i] = u[i]*h[i]
127            elif ( x[i] <= x3 ):
128                #print 'here x3=', x3
129                u[i] = u3
130                h[i] = h3
131                uh[i] = u[i]*h[i]
132            elif ( x[i] < x2 ):
133                #print 'here x2=', x2
134                u[i] = u2
135                h[i] = h2
136                uh[i] = u[i]*h[i]
137            else:
138                #print 'here the last section'
139                u[i] = 0.0 
140                h[i] = h0
141                uh[i] = u[i]*h[i]
142
143        return h , uh, u
Note: See TracBrowser for help on using the repository browser.