source: anuga_core/source/anuga/examples/run_sw_rectangle.py @ 4026

Last change on this file since 4026 was 4026, checked in by jack, 17 years ago

Moved the vpython visualiser to obsolete_code and cleared out things that call it from other code.

File size: 5.2 KB
Line 
1#!/usr/bin/env python
2"""
3  Example showing improved wet dry limiting.
4  This script runs the same simulation twice and stores the outputs in
5  old_limiter.sww and new_limiter.sww respectively.
6 
7  Authors: Steve Roberts
8  October
9"""
10
11#-----------------------------------------------------------------
12# Common structures
13#-----------------------------------------------------------------
14import time
15from Numeric import array
16from anuga.shallow_water import Domain
17from anuga.shallow_water import Reflective_boundary
18from anuga.abstract_2d_finite_volumes.mesh_factory import rectangular
19
20class Set_Stage:
21    """Set an initial condition with constant water height, for x<x0
22    """
23
24    def __init__(self, x0=0.25, x1=0.75, y0=0.0, y1=1.0, h=5.0, h0=0.0):
25        self.x0 = x0
26        self.x1 = x1
27        self.y0 = y0
28        self.y1 = y1
29        self.= h
30        self.h0 = h0
31
32    def __call__(self, x, y):
33        return self.h0 + self.h*((x>self.x0)&(x<self.x1)&(y>self.y0)&(y<self.y1))
34
35M = 30
36points, vertices, boundary = rectangular(M, M, len1 = 1.0, len2 = 1.0)
37
38yieldstep = 0.002
39finaltime = 0.06
40rect = [0.0, 0.0, 1.0, 1.0]
41
42
43#-----------------------------------------------------------------
44# Create domain for "old limiter" scenario
45#-----------------------------------------------------------------
46domain = Domain(points, vertices, boundary)
47domain.set_name('old_limiter_second_order')
48print 'Number of triangles =', len(domain)
49
50#-----------------------------------------------------------------
51# Boundaries and Initial conditions
52#-----------------------------------------------------------------
53R = Reflective_boundary(domain)
54domain.set_boundary( {'left': R, 'right': R, 'bottom': R, 'top': R} )
55domain.set_quantity('stage', Set_Stage(0.2, 0.4, 0.25, 0.75, 1.0, 0.00))
56
57#-----------------------------------------------------------------
58# Values for old limiter
59#-----------------------------------------------------------------
60domain.default_order = 2
61domain.beta_w      = 0.9
62domain.beta_w_dry  = 0.9
63domain.beta_uh     = 0.9
64domain.beta_uh_dry = 0.9
65domain.beta_vh     = 0.9
66domain.beta_vh_dry = 0.9
67
68#-----------------------------------------------------------------
69# Evolve
70#-----------------------------------------------------------------
71t0 = time.time()
72for t in domain.evolve(yieldstep = yieldstep, finaltime = finaltime):
73    domain.write_time()
74
75print 'That took %.2f seconds' %(time.time()-t0)
76print 'Note the small timesteps and the irregular flow'
77#raw_input('press return to continue')
78
79
80#-----------------------------------------------------------------
81# Create domain for "new limiter" scenario (2 order)
82#-----------------------------------------------------------------
83domain = Domain(points, vertices, boundary)
84domain.set_name('new_limiter_second_order')
85print 'Number of triangles =', len(domain)
86
87#-----------------------------------------------------------------
88# Boundaries and Initial conditions
89#-----------------------------------------------------------------
90R = Reflective_boundary(domain)
91domain.set_boundary( {'left': R, 'right': R, 'bottom': R, 'top': R} )
92domain.set_quantity('stage', Set_Stage(0.2, 0.4, 0.25, 0.75, 1.0, 0.00))
93
94#-----------------------------------------------------------------
95# Values for new limiter
96#-----------------------------------------------------------------
97domain.set_default_order(2)
98domain.minimum_allowed_height = 0.001
99domain.beta_w      = 1.0
100domain.beta_w_dry  = 0.2
101domain.beta_uh     = 1.0
102domain.beta_uh_dry = 0.2
103domain.beta_vh     = 1.0
104domain.beta_vh_dry = 0.2
105
106#-----------------------------------------------------------------
107# Evolve
108#-----------------------------------------------------------------
109t0 = time.time()
110for t in domain.evolve(yieldstep = yieldstep, finaltime = finaltime):
111    domain.write_time()
112
113print 'That took %.2f seconds' %(time.time()-t0)
114print 'Note more uniform and large timesteps'
115#raw_input('press return to continue')
116
117
118#-----------------------------------------------------------------
119# Create domain for "new limiter" scenario (1 order)
120#-----------------------------------------------------------------
121domain = Domain(points, vertices, boundary)
122domain.set_name('new_limiter_first_order')
123print 'Number of triangles =', len(domain)
124
125#-----------------------------------------------------------------
126# Boundaries and Initial conditions
127#-----------------------------------------------------------------
128R = Reflective_boundary(domain)
129domain.set_boundary( {'left': R, 'right': R, 'bottom': R, 'top': R} )
130domain.set_quantity('stage', Set_Stage(0.2, 0.4, 0.25, 0.75, 1.0, 0.00))
131
132#-----------------------------------------------------------------
133# Values for new limiter first order
134#-----------------------------------------------------------------
135domain.set_default_order(1)
136domain.minimum_allowed_height = 0.001
137domain.beta_w      = 1.0
138domain.beta_w_dry  = 0.2
139domain.beta_uh     = 1.0
140domain.beta_uh_dry = 0.2
141domain.beta_vh     = 1.0
142domain.beta_vh_dry = 0.2
143
144#-----------------------------------------------------------------
145# Evolve
146#-----------------------------------------------------------------
147t0 = time.time()
148for t in domain.evolve(yieldstep = yieldstep, finaltime = finaltime):
149    domain.write_time()
150
151print 'That took %.2f seconds' %(time.time()-t0)
152
153
Note: See TracBrowser for help on using the repository browser.