source: anuga_work/development/anuga_1d/dry_dam_vel.py @ 7818

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

Added Padarn's (2008/2009 summer student) files

File size: 2.1 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       
81domain = Domain(points)
82   
83domain.set_quantity('stage', stage)
84domain.set_boundary({'exterior': Reflective_boundary(domain)})
85domain.order = 2
86domain.set_timestepping_method('euler')
87domain.set_CFL(1.0)
88domain.set_limiter("vanleer")
89#domain.h0=0.0001
90
91t0 = time.time()
92
93for t in domain.evolve(yieldstep = yieldstep, finaltime = finaltime):
94    domain.write_time()
95
96print 'end'
97
98
Note: See TracBrowser for help on using the repository browser.