1 | """Example of shallow water wave equation. |
---|
2 | |
---|
3 | This is called Netherlands because it shows a dam with a gap in it and |
---|
4 | stylised housed behind it and below the water surface. |
---|
5 | |
---|
6 | """ |
---|
7 | |
---|
8 | ###################### |
---|
9 | # Module imports |
---|
10 | # |
---|
11 | from shallow_water import Domain, Reflective_boundary, Dirichlet_boundary,\ |
---|
12 | Transmissive_boundary, Time_boundary, Constant_height |
---|
13 | |
---|
14 | from mesh_factory import from_polyfile, rectangular |
---|
15 | from Numeric import array |
---|
16 | from math import sqrt |
---|
17 | #from least_squares import Interpolation |
---|
18 | from anuga.fit_interpolate.interpolate import Interpolate |
---|
19 | |
---|
20 | |
---|
21 | print 'Creating domain' |
---|
22 | #data_points, _, data_values = from_polyfile('cornell_room_medres') |
---|
23 | #points, triangles, values = from_polyfile('hires2') |
---|
24 | data_points, _, data_values = from_polyfile('hires2') |
---|
25 | |
---|
26 | |
---|
27 | #Regrid onto numerically stable mesh |
---|
28 | # |
---|
29 | #Compute regular mesh based on resolution and extent of data |
---|
30 | data_points = array(data_points) |
---|
31 | pmax = max(data_points) |
---|
32 | pmin = min(data_points) |
---|
33 | |
---|
34 | M = len(data_points) |
---|
35 | |
---|
36 | N = int(0.8*sqrt(M)) |
---|
37 | |
---|
38 | #print N |
---|
39 | |
---|
40 | mesh_points, vertices, boundary = rectangular(N, N, |
---|
41 | len1=pmax[0]-pmin[0], |
---|
42 | len2=pmax[1]-pmin[1], |
---|
43 | origin = pmin) |
---|
44 | |
---|
45 | |
---|
46 | #Compute smooth surface on new mesh based on values from old (regrid) |
---|
47 | print 'Interp' |
---|
48 | interp = Interpolate(mesh_points, vertices, alpha=10) |
---|
49 | mesh_values = interp.fit( data_points, data_values) # this has not been tested |
---|
50 | print 'Len mesh values', len(mesh_values) |
---|
51 | print 'Len mesh points', len(mesh_points) |
---|
52 | |
---|
53 | |
---|
54 | #Create shallow water domain |
---|
55 | print 'Creating domain' |
---|
56 | domain = Domain(mesh_points, vertices) #, boundary) |
---|
57 | |
---|
58 | domain.check_integrity() |
---|
59 | domain.default_order = 2 |
---|
60 | domain.smooth = True |
---|
61 | domain.reduction = min #Looks a lot better on top of steep slopes |
---|
62 | |
---|
63 | print "Number of triangles = ", len(domain) |
---|
64 | |
---|
65 | domain.store = True #Store for visualisation purposes |
---|
66 | domain.format = 'sww' #Native netcdf visualisation format |
---|
67 | import sys, os |
---|
68 | root, ext = os.path.splitext(sys.argv[0]) |
---|
69 | if domain.smooth is True: |
---|
70 | s = 'smooth' |
---|
71 | else: |
---|
72 | s = 'nonsmooth' |
---|
73 | domain.set_name(root + '_' + s) |
---|
74 | |
---|
75 | #Set bed-slope and friction |
---|
76 | manning = 0.0 |
---|
77 | |
---|
78 | print 'Field values' |
---|
79 | domain.set_quantity('elevation', mesh_values) |
---|
80 | domain.set_quantity('friction', manning) |
---|
81 | |
---|
82 | |
---|
83 | ###################### |
---|
84 | # Boundary conditions |
---|
85 | # |
---|
86 | print 'Boundaries' |
---|
87 | Br = Reflective_boundary(domain) |
---|
88 | domain.set_boundary({'exterior': Br}) |
---|
89 | |
---|
90 | |
---|
91 | |
---|
92 | ###################### |
---|
93 | #Initial condition |
---|
94 | # |
---|
95 | print 'Initial condition' |
---|
96 | |
---|
97 | #Define water height as a lump in one corner |
---|
98 | def height(x, y): |
---|
99 | from Numeric import zeros, Float |
---|
100 | |
---|
101 | N = len(x) |
---|
102 | assert N == len(y) |
---|
103 | |
---|
104 | xmin = min(x); xmax = max(x) |
---|
105 | ymin = min(y); ymax = max(y) |
---|
106 | |
---|
107 | xrange = xmax - xmin |
---|
108 | yrange = ymax - ymin |
---|
109 | |
---|
110 | z = zeros(N, Float) |
---|
111 | for i in range(N): |
---|
112 | if x[i] <= xmin + 0.25*xrange and y[i] <= ymin + 0.25*yrange: |
---|
113 | z[i] = 300 |
---|
114 | |
---|
115 | return z |
---|
116 | |
---|
117 | domain.set_quantity('stage', height) |
---|
118 | |
---|
119 | E = domain.quantities['elevation'].vertex_values |
---|
120 | L = domain.quantities['stage'].vertex_values |
---|
121 | domain.set_quantity('stage', E+L) |
---|
122 | |
---|
123 | #Evolve |
---|
124 | for t in domain.evolve(yieldstep = 0.05, finaltime = 5.0): |
---|
125 | domain.write_time() |
---|
126 | |
---|
127 | |
---|
128 | |
---|