1 | import os |
---|

2 | from math import sqrt, pi |
---|

3 | from shallow_water_vel_domain import * |
---|

4 | from Numeric import allclose, array, zeros, ones, Float, take, sqrt |
---|

5 | from config import g, epsilon |
---|

6 | |
---|

7 | |
---|

8 | h1 = 10.0 |
---|

9 | h0 = 0.0 |
---|

10 | |
---|

11 | def analytical_sol(C,t): |
---|

12 | |
---|

13 | #t = 0.0 # time (s) |
---|

14 | # gravity (m/s^2) |
---|

15 | #h1 = 10.0 # depth upstream (m) |
---|

16 | #h0 = 0.0 # depth downstream (m) |
---|

17 | L = 2000.0 # length of stream/domain (m) |
---|

18 | n = len(C) # number of cells |
---|

19 | |
---|

20 | u = zeros(n,Float) |
---|

21 | h = zeros(n,Float) |
---|

22 | x = C-3*L/4.0 |
---|

23 | |
---|

24 | |
---|

25 | for i in range(n): |
---|

26 | # Calculate Analytical Solution at time t > 0 |
---|

27 | u3 = 2.0/3.0*(sqrt(g*h1)+x[i]/t) |
---|

28 | h3 = 4.0/(9.0*g)*(sqrt(g*h1)-x[i]/(2.0*t))*(sqrt(g*h1)-x[i]/(2.0*t)) |
---|

29 | u3_ = 2.0/3.0*((x[i]+L/2.0)/t-sqrt(g*h1)) |
---|

30 | 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)) |
---|

31 | |
---|

32 | if ( x[i] <= -1*L/2.0+2*(-sqrt(g*h1)*t)): |
---|

33 | u[i] = 0.0 |
---|

34 | h[i] = h0 |
---|

35 | elif ( x[i] <= -1*L/2.0-(-sqrt(g*h1)*t)): |
---|

36 | u[i] = u3_ |
---|

37 | h[i] = h3_ |
---|

38 | |
---|

39 | elif ( x[i] <= -t*sqrt(g*h1) ): |
---|

40 | u[i] = 0.0 |
---|

41 | h[i] = h1 |
---|

42 | elif ( x[i] <= 2.0*t*sqrt(g*h1) ): |
---|

43 | u[i] = u3 |
---|

44 | h[i] = h3 |
---|

45 | else: |
---|

46 | u[i] = 0.0 |
---|

47 | h[i] = h0 |
---|

48 | |
---|

49 | return h , u*h, u |
---|

50 | |
---|

51 | print "TEST 1D-SOLUTION III -- DRY BED" |
---|

52 | |
---|

53 | def stage(x): |
---|

54 | y = zeros(len(x),Float) |
---|

55 | for i in range(len(x)): |
---|

56 | if x[i]<=L/4.0: |
---|

57 | y[i] = h0 |
---|

58 | elif x[i]<=3*L/4.0: |
---|

59 | y[i] = h1 |
---|

60 | else: |
---|

61 | y[i] = h0 |
---|

62 | return y |
---|

63 | |
---|

64 | |
---|

65 | import time |
---|

66 | |
---|

67 | finaltime = 2.0 |
---|

68 | yieldstep = 0.1 |
---|

69 | L = 2000.0 # Length of channel (m) |
---|

70 | |
---|

71 | k = 0 |
---|

72 | |
---|

73 | N = 800 |
---|

74 | print "Evaluating domain with %d cells" %N |
---|

75 | cell_len = L/N # Origin = 0.0 |
---|

76 | points = zeros(N+1,Float) |
---|

77 | |
---|

78 | for j in range(N+1): |
---|

79 | points[j] = j*cell_len |
---|

80 | |
---|

81 | boundary = {} |
---|

82 | boundary[0,0] = 'left' |
---|

83 | boundary[N-1,1] = 'right' |
---|

84 | |
---|

85 | domain = Domain(points, boundary = boundary) |
---|

86 | |
---|

87 | domain.set_quantity('stage', stage) |
---|

88 | |
---|

89 | Br = Reflective_boundary(domain) |
---|

90 | domain.set_boundary({'left': Br, 'right': Br}) |
---|

91 | domain.order = 2 |
---|

92 | domain.set_timestepping_method('euler') |
---|

93 | domain.set_CFL(1.0) |
---|

94 | domain.set_limiter("vanleer") |
---|

95 | #domain.h0=0.0001 |
---|

96 | |
---|

97 | t0 = time.time() |
---|

98 | |
---|

99 | for t in domain.evolve(yieldstep = yieldstep, finaltime = finaltime): |
---|

100 | domain.write_time() |
---|

101 | |
---|

102 | print 'end' |
---|

103 | |
---|

104 | |
---|