1 | """Validation study of Merimbula lake using Pyvolution. |
---|
2 | Example of shallow water wave equation applied to |
---|
3 | Malpasset dam break simulation. |
---|
4 | |
---|
5 | Copyright 2004 |
---|
6 | Christopher Zoppou, Stephen Roberts |
---|
7 | Australian National University |
---|
8 | |
---|
9 | Specific methods pertaining to the 2D shallow water equation |
---|
10 | are imported from shallow_water |
---|
11 | for use with the generic finite volume framework |
---|
12 | |
---|
13 | Conserved quantities are h, uh and vh stored as elements 0, 1 and 2 in the |
---|
14 | numerical vector named conserved_quantities. |
---|
15 | |
---|
16 | Existence of file 'Hobart_mesh.tsh' is assumed. |
---|
17 | """ |
---|
18 | |
---|
19 | ############################### |
---|
20 | # Setup Path and import modules |
---|
21 | import sys |
---|
22 | from os import sep, path |
---|
23 | sys.path.append('..'+sep+'pyvolution') |
---|
24 | |
---|
25 | from shallow_water import Domain, Reflective_boundary, File_boundary,\ |
---|
26 | Dirichlet_boundary, Transmissive_boundary, Constant_height |
---|
27 | from pmesh2domain import pmesh_to_domain_instance |
---|
28 | from util import inside_polygon, read_polygon |
---|
29 | |
---|
30 | ###################### |
---|
31 | # Domain |
---|
32 | filename = 'Hobart_mesh.tsh' |
---|
33 | yieldstep = 1 |
---|
34 | finaltime = 100 |
---|
35 | |
---|
36 | print 'Creating domain from', filename |
---|
37 | domain = pmesh_to_domain_instance(filename, Domain) |
---|
38 | print 'Number of triangles = ', len(domain) |
---|
39 | print 'Extent = ', domain.get_extent() |
---|
40 | |
---|
41 | |
---|
42 | domain.default_order = 1 |
---|
43 | domain.smooth = True |
---|
44 | |
---|
45 | #domain.visualise = True |
---|
46 | |
---|
47 | |
---|
48 | #------------------------------ |
---|
49 | # Boundary Conditions |
---|
50 | tags = {} |
---|
51 | tags['external'] = Reflective_boundary(domain) |
---|
52 | domain.set_boundary(tags) |
---|
53 | |
---|
54 | #----------------- |
---|
55 | #Initial condition |
---|
56 | # Throughout |
---|
57 | domain.set_quantity('stage', 0.) |
---|
58 | |
---|
59 | # Within polygon region |
---|
60 | # Set elevation to bed plus depth |
---|
61 | depth = 10. |
---|
62 | # Read polygon boundary |
---|
63 | p0 = read_polygon('Hobart_recent_source_zone.xya') |
---|
64 | # Get all centroid coordinates |
---|
65 | X = domain.get_centroid_coordinates() |
---|
66 | # Find triangle indices which are within polygon boundary only |
---|
67 | indices = inside_polygon(X, p0) |
---|
68 | # Get the bed elevation at the centroid of every triangle in polygon region |
---|
69 | |
---|
70 | zp = domain.get_quantity('elevation', location='centroids', indexes=indices) |
---|
71 | # Set the stage to bed elevation plus depth for all triangle vertices in the polygon region |
---|
72 | |
---|
73 | domain.set_quantity('stage', zp+depth, location='centroids', indexes=indices) |
---|
74 | |
---|
75 | |
---|
76 | #Alternative way which also works (I used it for testing - OMN) |
---|
77 | # |
---|
78 | |
---|
79 | #z = domain.get_quantity('elevation') |
---|
80 | #from copy import copy |
---|
81 | #w = copy(z) |
---|
82 | #for i in indices: |
---|
83 | # w[i, :] = z[i, :] + depth |
---|
84 | #domain.set_quantity('stage', w) |
---|
85 | |
---|
86 | |
---|
87 | #------------------------------------- |
---|
88 | # Provide file name for storing output |
---|
89 | domain.store = True |
---|
90 | domain.format = 'sww' |
---|
91 | domain.filename = 'Hobart_first_order' |
---|
92 | |
---|
93 | #--------------------------------------------------------- |
---|
94 | #Decide which quantities are to be stored at each timestep |
---|
95 | domain.quantities_to_be_stored = ['stage', 'xmomentum', 'ymomentum'] |
---|
96 | |
---|
97 | #---------------------------- |
---|
98 | # Friction |
---|
99 | domain.set_quantity('friction', 0.05) |
---|
100 | |
---|
101 | ###################### |
---|
102 | #Evolution |
---|
103 | import time |
---|
104 | t0 = time.time() |
---|
105 | for t in domain.evolve(yieldstep = yieldstep, finaltime = finaltime): |
---|
106 | domain.write_time() |
---|
107 | |
---|
108 | print 'That took %.2f seconds' %(time.time()-t0) |
---|
109 | |
---|
110 | |
---|