1 | """Validation of the AnuGA implementation of the shallow water wave equation. |
---|
2 | |
---|
3 | This script sets up LWRU2 benchmark with initial condition stated |
---|
4 | |
---|
5 | See also |
---|
6 | |
---|
7 | http://www.cee.cornell.edu/longwave/index.cfm?page=benchmark&problem=2 |
---|
8 | |
---|
9 | Depth at western boundary is d = 13.5 cm |
---|
10 | """ |
---|
11 | |
---|
12 | # Module imports |
---|
13 | from Numeric import array, zeros, Float, allclose |
---|
14 | |
---|
15 | from anuga.shallow_water import Domain |
---|
16 | from anuga.shallow_water import Reflective_boundary |
---|
17 | from anuga.shallow_water import Transmissive_Momentum_Set_Stage_boundary |
---|
18 | from anuga.abstract_2d_finite_volumes.util import file_function |
---|
19 | |
---|
20 | import project |
---|
21 | |
---|
22 | |
---|
23 | #------------------------- |
---|
24 | # Create Domain |
---|
25 | #------------------------- |
---|
26 | |
---|
27 | use_variable_mesh = True #Use large variable mesh generated by create_mesh.py |
---|
28 | #use_variable_mesh = False #Use small structured mesh for speed |
---|
29 | |
---|
30 | if use_variable_mesh is True: |
---|
31 | print 'Creating domain from', project.mesh_filename |
---|
32 | |
---|
33 | domain = Domain(project.mesh_filename, use_cache=True, verbose=True) |
---|
34 | else: |
---|
35 | print 'Creating regular mesh' |
---|
36 | N = 150 |
---|
37 | points, vertices, boundary = rectangular_cross(N, N/5*3, |
---|
38 | len1=5.448, len2=3.402) |
---|
39 | print 'Creating domain' |
---|
40 | |
---|
41 | #domain = Domain(points, vertices, boundary) |
---|
42 | domain = cache(Domain, (points, vertices, boundary)) |
---|
43 | |
---|
44 | |
---|
45 | domain.set_name('lwru2') |
---|
46 | domain.set_default_order(2) |
---|
47 | domain.set_minimum_storable_height(0.001) |
---|
48 | |
---|
49 | #domain.check_integrity() |
---|
50 | print 'Number of triangles = ', len(domain) |
---|
51 | print 'The extent is ', domain.get_extent() |
---|
52 | print domain.statistics() |
---|
53 | |
---|
54 | |
---|
55 | #------------------------- |
---|
56 | # Initial Conditions |
---|
57 | #------------------------- |
---|
58 | domain.set_quantity('friction', 0.0) |
---|
59 | domain.set_quantity('stage', 0.0) |
---|
60 | domain.set_quantity('elevation', |
---|
61 | filename = project.bathymetry_filename[:-4] + '.pts', |
---|
62 | alpha = 0.02, |
---|
63 | verbose = True, |
---|
64 | use_cache = True) |
---|
65 | |
---|
66 | #------------------------- |
---|
67 | # Boundary Conditions |
---|
68 | #------------------------- |
---|
69 | Br = Reflective_boundary(domain) |
---|
70 | |
---|
71 | function = file_function(project.boundary_filename[:-4] + '.tms', |
---|
72 | domain, |
---|
73 | verbose = True) |
---|
74 | |
---|
75 | Bts = Transmissive_Momentum_Set_Stage_boundary(domain, function) |
---|
76 | |
---|
77 | #Set boundary conditions |
---|
78 | if use_variable_mesh is True: |
---|
79 | domain.set_boundary({'wave': Bts, 'wall': Br}) |
---|
80 | else: |
---|
81 | domain.set_boundary({'left': Bts, 'right': Br, 'bottom': Br, 'top': Br}) |
---|
82 | |
---|
83 | |
---|
84 | #------------------------- |
---|
85 | # Evolve through time |
---|
86 | #------------------------- |
---|
87 | import time |
---|
88 | t0 = time.time() |
---|
89 | |
---|
90 | for t in domain.evolve(yieldstep = 0.05, finaltime = 22.5): |
---|
91 | domain.write_time() |
---|
92 | |
---|
93 | print 'That took %.2f seconds' %(time.time()-t0) |
---|