1 | """Validation of the AnuGA implementation of the shallow water wave equation. |
---|
2 | |
---|
3 | This script sets up 2D circular island..... |
---|
4 | |
---|
5 | """ |
---|
6 | |
---|
7 | # Module imports |
---|
8 | from anuga.shallow_water import Domain |
---|
9 | from anuga.shallow_water import Reflective_boundary |
---|
10 | from anuga.shallow_water import Transmissive_Momentum_Set_Stage_boundary |
---|
11 | from anuga.shallow_water import Dirichlet_boundary, Time_boundary |
---|
12 | from anuga.abstract_2d_finite_volumes.util import file_function |
---|
13 | from anuga.pmesh.mesh_interface import create_mesh_from_regions |
---|
14 | #from cmath import cos,sin,cosh,pi |
---|
15 | from math import cos,pi,sin,tan,sqrt |
---|
16 | from Numeric import array, zeros, Float, allclose |
---|
17 | |
---|
18 | #import project |
---|
19 | |
---|
20 | #------------------------- |
---|
21 | # Create Domain from mesh |
---|
22 | #------------------------- |
---|
23 | |
---|
24 | |
---|
25 | def prepare_timeboundary(filename): |
---|
26 | """Converting benchmark 2 time series to NetCDF tms file. |
---|
27 | This is a 'throw-away' code taylor made for files like |
---|
28 | 'Benchmark_2_input.txt' from the LWRU2 benchmark |
---|
29 | """ |
---|
30 | |
---|
31 | textversion = filename[:-4] + '.txt' |
---|
32 | |
---|
33 | print 'Preparing time boundary from %s' %textversion |
---|
34 | from Numeric import array |
---|
35 | |
---|
36 | fid = open(textversion) |
---|
37 | |
---|
38 | #Skip first line |
---|
39 | line = fid.readline() |
---|
40 | |
---|
41 | #Read remaining lines |
---|
42 | lines = fid.readlines() |
---|
43 | fid.close() |
---|
44 | |
---|
45 | |
---|
46 | N = len(lines) |
---|
47 | T = zeros(N, Float) #Time |
---|
48 | Q = zeros(N, Float) #Values |
---|
49 | |
---|
50 | for i, line in enumerate(lines): |
---|
51 | fields = line.split() |
---|
52 | print 'fields',i,line |
---|
53 | |
---|
54 | T[i] = float(fields[0]) |
---|
55 | Q[i] = float(fields[1]) |
---|
56 | |
---|
57 | |
---|
58 | #Create tms file |
---|
59 | from Scientific.IO.NetCDF import NetCDFFile |
---|
60 | |
---|
61 | print 'Writing to', filename |
---|
62 | fid = NetCDFFile(filename[:-4] + '.tms', 'w') |
---|
63 | |
---|
64 | fid.institution = 'Geoscience Australia' |
---|
65 | fid.description = 'Input wave for Benchmark 2' |
---|
66 | fid.starttime = 0.0 |
---|
67 | fid.createDimension('number_of_timesteps', len(T)) |
---|
68 | fid.createVariable('time', Float, ('number_of_timesteps',)) |
---|
69 | fid.variables['time'][:] = T |
---|
70 | |
---|
71 | fid.createVariable('stage', Float, ('number_of_timesteps',)) |
---|
72 | fid.variables['stage'][:] = Q[:] |
---|
73 | |
---|
74 | fid.createVariable('xmomentum', Float, ('number_of_timesteps',)) |
---|
75 | fid.variables['xmomentum'][:] = 0.0 |
---|
76 | |
---|
77 | fid.createVariable('ymomentum', Float, ('number_of_timesteps',)) |
---|
78 | fid.variables['ymomentum'][:] = 0.0 |
---|
79 | |
---|
80 | fid.close() |
---|
81 | |
---|
82 | |
---|
83 | |
---|
84 | |
---|
85 | poly_all= [[0,0],[0,25],[30,25],[30,0]] |
---|
86 | poly=[[9,10],[9,17.6],[16.6,17.6],[16.6,10]] |
---|
87 | internal_poly=[[poly,1]] |
---|
88 | create_mesh_from_regions(poly_all, |
---|
89 | boundary_tags={'wall': [0,1,2],'wave': [3]}, |
---|
90 | maximum_triangle_area=10, |
---|
91 | interior_regions=internal_poly, |
---|
92 | filename='temp.msh', |
---|
93 | use_cache=False, |
---|
94 | verbose=True) |
---|
95 | |
---|
96 | domain = Domain('temp.msh', use_cache=True, verbose=True) |
---|
97 | print domain.statistics() |
---|
98 | |
---|
99 | #------------------------- |
---|
100 | # Initial Conditions |
---|
101 | #------------------------- |
---|
102 | domain.set_quantity('friction', 0.0) |
---|
103 | water_height = 0.32 |
---|
104 | domain.set_quantity('stage', water_height) |
---|
105 | |
---|
106 | def test(): |
---|
107 | fileName = "test.csv" |
---|
108 | file = open(fileName,"w") |
---|
109 | file.write("x,y,elevation \n") |
---|
110 | x = 12.96 |
---|
111 | y = 13.80 |
---|
112 | for i in range(0,30): |
---|
113 | for j in range(0,25): |
---|
114 | a = (x-i)**2 |
---|
115 | b = (y-j)**2 |
---|
116 | z = -sqrt((x-i)**2+(y-j)**2) |
---|
117 | z = z*.25+(0.8975) |
---|
118 | if z <0: |
---|
119 | z=0 |
---|
120 | if z > .625: |
---|
121 | z=0.625 |
---|
122 | |
---|
123 | print 'x,y,f',i,j,a,b,z |
---|
124 | file.write("%s, %s, %s\n" %(i, j, z)) |
---|
125 | file.close() |
---|
126 | |
---|
127 | test() |
---|
128 | |
---|
129 | domain.set_quantity('elevation',filename = "test.csv", alpha=0.1) |
---|
130 | |
---|
131 | |
---|
132 | #------------------------- |
---|
133 | # Set simulation parameters |
---|
134 | #------------------------- |
---|
135 | domain.set_name('test') # Name of output sww file |
---|
136 | domain.set_default_order(1) # Apply second order scheme |
---|
137 | domain.set_all_limiters(0.9) # Max second order scheme (old lim) |
---|
138 | domain.set_minimum_storable_height(0.01) # Don't store w < 0.001m |
---|
139 | domain.set_maximum_allowed_speed(0.1) # Allow a little runoff (0.1 is OK) |
---|
140 | domain.tight_slope_limiters = 1 |
---|
141 | domain.beta_h = 0.0 |
---|
142 | |
---|
143 | #Timings on AMD64-242 (beta_h=0) |
---|
144 | # tight_slope_limiters = 0: |
---|
145 | # 3035s - 3110s |
---|
146 | # tight_slope_limiters = 1: |
---|
147 | # 3000s - 3008s |
---|
148 | # |
---|
149 | # beta_h>0: In the order of 3200s |
---|
150 | |
---|
151 | #------------------------- |
---|
152 | # Boundary Conditions |
---|
153 | #------------------------- |
---|
154 | |
---|
155 | # Create boundary function from timeseries provided in file |
---|
156 | #function = file_function(project.boundary_filename, |
---|
157 | # domain, verbose=True) |
---|
158 | def waveform(t): |
---|
159 | return (sin(t*2*pi)) |
---|
160 | |
---|
161 | # Create and assign boundary objects |
---|
162 | boundary_filename="ts2a_edit_input_32" |
---|
163 | prepare_timeboundary(boundary_filename+'.txt') |
---|
164 | |
---|
165 | function = file_function(boundary_filename+'.tms', |
---|
166 | domain, verbose=True) |
---|
167 | |
---|
168 | # Create and assign boundary objects |
---|
169 | Bts = Transmissive_Momentum_Set_Stage_boundary(domain, function) |
---|
170 | |
---|
171 | #Bts = Transmissive_Momentum_Set_Stage_boundary(domain, waveform) |
---|
172 | Bw = Time_boundary(domain=domain, |
---|
173 | f=lambda t: [(0.1*sin(t*2*pi)), 0.0, 0.0]) |
---|
174 | Br = Reflective_boundary(domain) |
---|
175 | Bd = Dirichlet_boundary([water_height,0,0]) |
---|
176 | domain.set_boundary({'wave': Bts, 'wall': Bd}) |
---|
177 | |
---|
178 | #------------------------- |
---|
179 | # Evolve through time |
---|
180 | #------------------------- |
---|
181 | import time |
---|
182 | t0 = time.time() |
---|
183 | |
---|
184 | for t in domain.evolve(yieldstep = 5, finaltime = 75): |
---|
185 | domain.write_time() |
---|
186 | # domain.write_time(track_speeds=False) |
---|
187 | |
---|
188 | print 'That took %.2f seconds' %(time.time()-t0) |
---|