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

Copied John Jakeman's and Sudi Mungkasi's 1d shallow water code

File size:
3.8 KB

Line  

1  

2  

3  class 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 = 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=zc0*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/c0u2/c02.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 = C3*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/(h2h0)*(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.0u2/z) 

78  x3=(u2c2)*t 

79  x2=z*t 

80  x1=c1*t 

81  x1_ = 1*L/2.0x1 

82  x2_ = 1*L/2.0+2*x1 

83  #x3_ = 1*L/2.0x3 

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)/tsqrt(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.