source: anuga_work/development/2010-projects/anuga_1d/channel/profile_channel.py @ 7855

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

Changing for loop to numpy.where

File size: 2.0 KB
Line 
1#! /usr/bin/python
2
3# To change this template, choose Tools | Templates
4# and open the template in the editor.
5
6__author__="steve"
7__date__ ="$15/06/2010 5:02:36 PM$"
8
9
10import os
11from math import sqrt, pow, pi
12
13import numpy
14
15from anuga_1d.channel.channel_domain import *
16from anuga_1d.config import g, epsilon
17from anuga_1d.base.generic_mesh import uniform_mesh
18
19
20def run_evolve():
21
22
23    print "Channel Flow 1 Padarn Test"
24
25    # Define functions for initial quantities
26    def initial_area(x):
27        return 1.4691*width(x)
28
29    def width(x):
30        x1=(x/1000)*(x/1000)
31        x2=x1*(x/1000)
32        x3=x2*(x/1000)
33        return 10-64*(x1-2*x2+x3)
34
35    def bed(x):
36        y = numpy.ones(len(x),numpy.float)
37
38        return numpy.where( (x<525) & (x>475),y,0.0)
39
40
41    def initial_discharge(x):
42        return 20
43
44    import time
45
46    # Set final time and yield time for simulation
47    finaltime = 50.0
48    yieldstep = 10.0
49
50    # Length of channel (m)
51    L = 1000.0
52    # Define the number of cells
53    N = 200
54
55    # Create domain with centroid points as defined above
56    domain = Domain(*uniform_mesh(N))
57
58
59    # Set initial values of quantities - default to zero
60    domain.set_quantity('area', initial_area)
61    domain.set_quantity('width',width)
62    domain.set_quantity('elevation',bed)
63    domain.set_quantity('discharge',initial_discharge)
64
65    # Set boundry type, order, timestepping method and limiter
66    Bd = Dirichlet_boundary([14,20,0,1.4,20/14,9,1.4])
67    domain.set_boundary({'left': Bd , 'right' : Bd })
68
69
70    domain.order = 2
71    domain.set_timestepping_method('rk2')
72    domain.set_CFL(1.0)
73    domain.set_limiter("vanleer")
74    #domain.set_limiter("minmod")
75    #domain.h0=0.0001
76
77    # Start timer
78    t0 = time.time()
79
80    for t in domain.evolve(yieldstep = yieldstep, finaltime = finaltime):
81        domain.write_time()
82
83
84
85import cProfile
86cProfile.run('run_evolve()', 'evolve_prof')
87
88
89import pstats
90p = pstats.Stats('evolve_prof')
91
92#p.strip_dirs().sort_stats(-1).print_stats(20)
93
94p.sort_stats('cumulative').print_stats(30)
95
96print p
97
98
99
100
Note: See TracBrowser for help on using the repository browser.