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 | #----------------------------------------------------------------- |
---|
14 | import time |
---|
15 | from Numeric import array |
---|
16 | from anuga.shallow_water import Domain |
---|
17 | from anuga.shallow_water import Reflective_boundary |
---|
18 | from anuga.abstract_2d_finite_volumes.mesh_factory import rectangular |
---|
19 | |
---|
20 | class 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 = 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 | |
---|
35 | M = 30 |
---|
36 | points, vertices, boundary = rectangular(M, M, len1 = 1.0, len2 = 1.0) |
---|
37 | |
---|
38 | yieldstep = 0.002 |
---|
39 | finaltime = 0.06 |
---|
40 | rect = [0.0, 0.0, 1.0, 1.0] |
---|
41 | |
---|
42 | |
---|
43 | #----------------------------------------------------------------- |
---|
44 | # Create domain for "old limiter" scenario |
---|
45 | #----------------------------------------------------------------- |
---|
46 | domain = Domain(points, vertices, boundary) |
---|
47 | domain.set_name('old_limiter_second_order') |
---|
48 | print 'Number of triangles =', len(domain) |
---|
49 | |
---|
50 | #----------------------------------------------------------------- |
---|
51 | # Boundaries and Initial conditions |
---|
52 | #----------------------------------------------------------------- |
---|
53 | R = Reflective_boundary(domain) |
---|
54 | domain.set_boundary( {'left': R, 'right': R, 'bottom': R, 'top': R} ) |
---|
55 | domain.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 | #----------------------------------------------------------------- |
---|
60 | domain.default_order = 2 |
---|
61 | domain.beta_w = 0.9 |
---|
62 | domain.beta_w_dry = 0.9 |
---|
63 | domain.beta_uh = 0.9 |
---|
64 | domain.beta_uh_dry = 0.9 |
---|
65 | domain.beta_vh = 0.9 |
---|
66 | domain.beta_vh_dry = 0.9 |
---|
67 | |
---|
68 | #----------------------------------------------------------------- |
---|
69 | # Evolve |
---|
70 | #----------------------------------------------------------------- |
---|
71 | t0 = time.time() |
---|
72 | for t in domain.evolve(yieldstep = yieldstep, finaltime = finaltime): |
---|
73 | domain.write_time() |
---|
74 | |
---|
75 | print 'That took %.2f seconds' %(time.time()-t0) |
---|
76 | print '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 | #----------------------------------------------------------------- |
---|
83 | domain = Domain(points, vertices, boundary) |
---|
84 | domain.set_name('new_limiter_second_order') |
---|
85 | print 'Number of triangles =', len(domain) |
---|
86 | |
---|
87 | #----------------------------------------------------------------- |
---|
88 | # Boundaries and Initial conditions |
---|
89 | #----------------------------------------------------------------- |
---|
90 | R = Reflective_boundary(domain) |
---|
91 | domain.set_boundary( {'left': R, 'right': R, 'bottom': R, 'top': R} ) |
---|
92 | domain.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 | #----------------------------------------------------------------- |
---|
97 | domain.set_default_order(2) |
---|
98 | domain.minimum_allowed_height = 0.001 |
---|
99 | domain.beta_w = 1.0 |
---|
100 | domain.beta_w_dry = 0.2 |
---|
101 | domain.beta_uh = 1.0 |
---|
102 | domain.beta_uh_dry = 0.2 |
---|
103 | domain.beta_vh = 1.0 |
---|
104 | domain.beta_vh_dry = 0.2 |
---|
105 | |
---|
106 | #----------------------------------------------------------------- |
---|
107 | # Evolve |
---|
108 | #----------------------------------------------------------------- |
---|
109 | t0 = time.time() |
---|
110 | for t in domain.evolve(yieldstep = yieldstep, finaltime = finaltime): |
---|
111 | domain.write_time() |
---|
112 | |
---|
113 | print 'That took %.2f seconds' %(time.time()-t0) |
---|
114 | print '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 | #----------------------------------------------------------------- |
---|
121 | domain = Domain(points, vertices, boundary) |
---|
122 | domain.set_name('new_limiter_first_order') |
---|
123 | print 'Number of triangles =', len(domain) |
---|
124 | |
---|
125 | #----------------------------------------------------------------- |
---|
126 | # Boundaries and Initial conditions |
---|
127 | #----------------------------------------------------------------- |
---|
128 | R = Reflective_boundary(domain) |
---|
129 | domain.set_boundary( {'left': R, 'right': R, 'bottom': R, 'top': R} ) |
---|
130 | domain.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 | #----------------------------------------------------------------- |
---|
135 | domain.set_default_order(1) |
---|
136 | domain.minimum_allowed_height = 0.001 |
---|
137 | domain.beta_w = 1.0 |
---|
138 | domain.beta_w_dry = 0.2 |
---|
139 | domain.beta_uh = 1.0 |
---|
140 | domain.beta_uh_dry = 0.2 |
---|
141 | domain.beta_vh = 1.0 |
---|
142 | domain.beta_vh_dry = 0.2 |
---|
143 | |
---|
144 | #----------------------------------------------------------------- |
---|
145 | # Evolve |
---|
146 | #----------------------------------------------------------------- |
---|
147 | t0 = time.time() |
---|
148 | for t in domain.evolve(yieldstep = yieldstep, finaltime = finaltime): |
---|
149 | domain.write_time() |
---|
150 | |
---|
151 | print 'That took %.2f seconds' %(time.time()-t0) |
---|
152 | |
---|
153 | |
---|