# source:anuga_work/development/anuga_1d/dry_dam_vel.py@6452

Last change on this file since 6452 was 6452, checked in by steve, 15 years ago

Just storing any changes

File size: 2.2 KB
Line
1import os
2from math import sqrt, pi
3from shallow_water_vel_domain import *
4from Numeric import allclose, array, zeros, ones, Float, take, sqrt
5from config import g, epsilon
6
7
8h1 = 10.0
9h0 = 0.0
10
11def 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
51print "TEST 1D-SOLUTION III -- DRY BED"
52
53def 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
65import time
66
67finaltime = 2.0
68yieldstep = 0.1
69L = 2000.0     # Length of channel (m)
70
71k = 0
72
73N = 800
74print "Evaluating domain with %d cells" %N
75cell_len = L/N # Origin = 0.0
76points = zeros(N+1,Float)
77
78for j in range(N+1):
79    points[j] = j*cell_len
80
81boundary = {}
82boundary[0,0] = 'left'
83boundary[N-1,1] = 'right'
84
85domain = Domain(points, boundary = boundary)
86
87domain.set_quantity('stage', stage)
88
89Br =  Reflective_boundary(domain)
90domain.set_boundary({'left':  Br, 'right': Br})
91domain.order = 2
92domain.set_timestepping_method('euler')
93domain.set_CFL(1.0)
94domain.set_limiter("vanleer")
95#domain.h0=0.0001
96
97t0 = time.time()
98
99for t in domain.evolve(yieldstep = yieldstep, finaltime = finaltime):
100    domain.write_time()
101
102print 'end'
103
104
Note: See TracBrowser for help on using the repository browser.