source: trunk/anuga_work/development/2010-projects/anuga_1d/sww/profile_dry_dam.py @ 7884

Last change on this file since 7884 was 7884, checked in by steve, 12 years ago

Moving 2010 project

File size: 2.7 KB
Line 
1import os
2from math import sqrt, pi
3import numpy
4import time
5#from Numeric import allclose, array, zeros, ones, Float, take, sqrt
6
7
8
9from anuga_1d.sww.sww_domain import *
10from anuga_1d.config import g, epsilon
11from anuga_1d.base.generic_mesh import uniform_mesh
12
13
14def run_evolve():
15
16    h1 = 10.0
17    h0 = 0.01
18
19    def analytical_sol(C,t):
20
21        #t  = 0.0     # time (s)
22        # gravity (m/s^2)
23        #h1 = 10.0    # depth upstream (m)
24        #h0 = 0.0     # depth downstream (m)
25        L = 2000.0   # length of stream/domain (m)
26        n = len(C)    # number of cells
27
28        u = zeros(n,Float)
29        h = zeros(n,Float)
30        x = C-3*L/4.0
31
32
33        for i in range(n):
34            # Calculate Analytical Solution at time t > 0
35            u3 = 2.0/3.0*(sqrt(g*h1)+x[i]/t)
36            h3 = 4.0/(9.0*g)*(sqrt(g*h1)-x[i]/(2.0*t))*(sqrt(g*h1)-x[i]/(2.0*t))
37            u3_ = 2.0/3.0*((x[i]+L/2.0)/t-sqrt(g*h1))
38            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))
39
40            if ( x[i] <= -1*L/2.0+2*(-sqrt(g*h1)*t)):
41                u[i] = 0.0
42                h[i] = h0
43            elif ( x[i] <= -1*L/2.0-(-sqrt(g*h1)*t)):
44                u[i] = u3_
45                h[i] = h3_
46
47            elif ( x[i] <= -t*sqrt(g*h1) ):
48                u[i] = 0.0
49                h[i] = h1
50            elif ( x[i] <= 2.0*t*sqrt(g*h1) ):
51                u[i] = u3
52                h[i] = h3
53            else:
54                u[i] = 0.0
55                h[i] = h0
56
57        return h , u*h, u
58
59
60
61    def stage(x):
62        import numpy
63
64        y = numpy.where( (x>=L/4.0) & (x<=3*L/4.0), h1 , h0)
65
66    #    for i in range(len(x)):
67    #        if x[i]<=L/4.0:
68    #            y[i] = h0
69    #        elif x[i]<=3*L/4.0:
70    #            y[i] = h1
71    #        else:
72    #            y[i] = h0
73        return y
74
75
76
77
78    print "TEST 1D-SOLUTION III -- DRY BED"
79
80
81    finaltime = 4.0
82    yieldstep = 0.1
83    L = 2000.0     # Length of channel (m)
84
85    k = 0
86
87    N = 800
88    print "Evaluating domain with %d cells" %N
89    domain = Domain(*uniform_mesh(N))
90
91    domain.set_quantity('stage', stage)
92
93    Br = Reflective_boundary(domain)
94
95    domain.set_boundary({'left': Br, 'right' : Br})
96    domain.order = 2
97    domain.set_timestepping_method('rk2')
98    domain.set_CFL(1.0)
99    domain.set_limiter("minmod")
100    #domain.h0=0.0001
101
102    t0 = time.time()
103
104    for t in domain.evolve(yieldstep = yieldstep, finaltime = finaltime):
105        domain.write_time()
106
107    print 'That took %.2f seconds' %(time.time()-t0)
108
109
110import cProfile
111cProfile.run('run_evolve()', 'dry_dam_prof')
112
113
114import pstats
115p = pstats.Stats('dry_dam_prof')
116
117#p.strip_dirs().sort_stats(-1).print_stats(20)
118
119p.sort_stats('cumulative').print_stats(30)
120
121
122
123
124
Note: See TracBrowser for help on using the repository browser.