1 | |
---|
2 | from anuga.shallow_water.shallow_water_domain import * |
---|
3 | from anuga.shallow_water.shallow_water_domain import Domain as Sww_domain |
---|
4 | |
---|
5 | |
---|
6 | ############################################################################## |
---|
7 | # Shallow Water Balanced Domain -- alternative implementation |
---|
8 | # |
---|
9 | # FIXME: Following the methods in CITE MODSIM PAPER |
---|
10 | # |
---|
11 | ############################################################################## |
---|
12 | |
---|
13 | class Domain(Sww_domain): |
---|
14 | |
---|
15 | def __init__(self, |
---|
16 | coordinates=None, |
---|
17 | vertices=None, |
---|
18 | boundary=None, |
---|
19 | tagged_elements=None, |
---|
20 | geo_reference=None, |
---|
21 | use_inscribed_circle=False, |
---|
22 | mesh_filename=None, |
---|
23 | use_cache=False, |
---|
24 | verbose=False, |
---|
25 | full_send_dict=None, |
---|
26 | ghost_recv_dict=None, |
---|
27 | starttime=0.0, |
---|
28 | processor=0, |
---|
29 | numproc=1, |
---|
30 | number_of_full_nodes=None, |
---|
31 | number_of_full_triangles=None): |
---|
32 | |
---|
33 | conserved_quantities = [ 'stage', 'xmomentum', 'ymomentum'] |
---|
34 | |
---|
35 | evolved_quantities = [ 'stage', 'xmomentum', 'ymomentum'] |
---|
36 | other_quantities = [ 'elevation', 'friction', 'height', |
---|
37 | 'xvelocity', 'yvelocity', 'x', 'y' ] |
---|
38 | |
---|
39 | |
---|
40 | Sww_domain.__init__(self, |
---|
41 | coordinates = coordinates, |
---|
42 | vertices = vertices, |
---|
43 | boundary = boundary, |
---|
44 | tagged_elements = tagged_elements, |
---|
45 | geo_reference = geo_reference, |
---|
46 | use_inscribed_circle = use_inscribed_circle, |
---|
47 | mesh_filename = mesh_filename, |
---|
48 | use_cache = use_cache, |
---|
49 | verbose = verbose, |
---|
50 | conserved_quantities = conserved_quantities, |
---|
51 | evolved_quantities = evolved_quantities, |
---|
52 | other_quantities = other_quantities, |
---|
53 | full_send_dict = full_send_dict, |
---|
54 | ghost_recv_dict = ghost_recv_dict, |
---|
55 | starttime = starttime, |
---|
56 | processor = processor, |
---|
57 | numproc = numproc, |
---|
58 | number_of_full_nodes = number_of_full_nodes, |
---|
59 | number_of_full_triangles = number_of_full_triangles) |
---|
60 | |
---|
61 | #--------------------- |
---|
62 | # set some defaults |
---|
63 | #--------------------- |
---|
64 | self.set_CFL(1.0) |
---|
65 | self.set_use_kinematic_viscosity(False) |
---|
66 | |
---|
67 | # Because gravity is treated within the flux function, |
---|
68 | # we remove it from the forcing terms. |
---|
69 | self.forcing_terms.remove(gravity) |
---|
70 | print 'Using shallow_water_balanced2 solver in /balanced_basic/' |
---|
71 | |
---|
72 | |
---|
73 | |
---|
74 | #----------------- |
---|
75 | # Flux computation |
---|
76 | #----------------- |
---|
77 | |
---|
78 | ## @brief Compute fluxes and timestep suitable for all volumes in domain. |
---|
79 | # @param domain The domain to calculate fluxes for. |
---|
80 | |
---|
81 | def compute_fluxes(domain): |
---|
82 | """Compute fluxes and timestep suitable for all volumes in domain. |
---|
83 | |
---|
84 | Compute total flux for each conserved quantity using "flux_function" |
---|
85 | |
---|
86 | Fluxes across each edge are scaled by edgelengths and summed up |
---|
87 | Resulting flux is then scaled by area and stored in |
---|
88 | explicit_update for each of the three conserved quantities |
---|
89 | stage, xmomentum and ymomentum |
---|
90 | |
---|
91 | The maximal allowable speed computed by the flux_function for each volume |
---|
92 | is converted to a timestep that must not be exceeded. The minimum of |
---|
93 | those is computed as the next overall timestep. |
---|
94 | |
---|
95 | Post conditions: |
---|
96 | domain.explicit_update is reset to computed flux values |
---|
97 | domain.timestep is set to the largest step satisfying all volumes. |
---|
98 | |
---|
99 | This wrapper calls the underlying C version of compute fluxes |
---|
100 | """ |
---|
101 | |
---|
102 | import sys |
---|
103 | from swb2_domain_ext import compute_fluxes_ext_central \ |
---|
104 | as compute_fluxes_ext |
---|
105 | |
---|
106 | # Shortcuts |
---|
107 | Stage = domain.quantities['stage'] |
---|
108 | Xmom = domain.quantities['xmomentum'] |
---|
109 | Ymom = domain.quantities['ymomentum'] |
---|
110 | Bed = domain.quantities['elevation'] |
---|
111 | |
---|
112 | timestep = float(sys.maxint) |
---|
113 | |
---|
114 | flux_timestep = compute_fluxes_ext(timestep, |
---|
115 | domain.epsilon, |
---|
116 | domain.H0, |
---|
117 | domain.g, |
---|
118 | domain.neighbours, |
---|
119 | domain.neighbour_edges, |
---|
120 | domain.normals, |
---|
121 | domain.edgelengths, |
---|
122 | domain.radii, |
---|
123 | domain.areas, |
---|
124 | domain.tri_full_flag, |
---|
125 | Stage.edge_values, |
---|
126 | Xmom.edge_values, |
---|
127 | Ymom.edge_values, |
---|
128 | Bed.edge_values, |
---|
129 | Stage.boundary_values, |
---|
130 | Xmom.boundary_values, |
---|
131 | Ymom.boundary_values, |
---|
132 | Stage.explicit_update, |
---|
133 | Xmom.explicit_update, |
---|
134 | Ymom.explicit_update, |
---|
135 | domain.already_computed_flux, |
---|
136 | domain.max_speed, |
---|
137 | int(domain.optimise_dry_cells), |
---|
138 | Stage.centroid_values, |
---|
139 | Bed.centroid_values, |
---|
140 | Bed.vertex_values) |
---|
141 | |
---|
142 | #import pdb |
---|
143 | #pdb.set_trace() |
---|
144 | |
---|
145 | domain.flux_timestep = flux_timestep |
---|
146 | |
---|
147 | |
---|
148 | def protect_against_infinitesimal_and_negative_heights(domain): |
---|
149 | """protect against infinitesimal heights and associated high velocities""" |
---|
150 | |
---|
151 | from swb2_domain_ext import protect |
---|
152 | #print'using swb2_protect_against ..' |
---|
153 | |
---|
154 | # shortcuts |
---|
155 | wc = domain.quantities['stage'].centroid_values |
---|
156 | zc = domain.quantities['elevation'].centroid_values |
---|
157 | zv = domain.quantities['elevation'].vertex_values |
---|
158 | xmomc = domain.quantities['xmomentum'].centroid_values |
---|
159 | ymomc = domain.quantities['ymomentum'].centroid_values |
---|
160 | |
---|
161 | protect(domain.minimum_allowed_height, domain.maximum_allowed_speed, |
---|
162 | domain.epsilon, wc, zc,zv, xmomc, ymomc) |
---|
163 | |
---|
164 | def conserved_values_to_evolved_values(self, q_cons, q_evol): |
---|
165 | """Mapping between conserved quantities and the evolved quantities. |
---|
166 | Used where we have a boundary condition which works with conserved |
---|
167 | quantities and we now want to use them for the new well balanced |
---|
168 | code using the evolved quantities |
---|
169 | |
---|
170 | Typically the old boundary condition will set the values in q_cons, |
---|
171 | |
---|
172 | q_evol on input will have the values of the evolved quantities at the |
---|
173 | edge point (useful as it provides values for evlevation). |
---|
174 | |
---|
175 | """ |
---|
176 | |
---|
177 | wc = q_cons[0] |
---|
178 | uhc = q_cons[1] |
---|
179 | vhc = q_cons[2] |
---|
180 | |
---|
181 | |
---|
182 | q_evol[0] = wc |
---|
183 | q_evol[1] = uhc |
---|
184 | q_evol[2] = vhc |
---|
185 | |
---|
186 | return q_evol |
---|
187 | |
---|
188 | def distribute_to_vertices_and_edges(self): |
---|
189 | """ Call correct module function """ |
---|
190 | |
---|
191 | #Shortcuts |
---|
192 | #W = self.quantities['stage'] |
---|
193 | #Z = self.quantities['elevation'] |
---|
194 | |
---|
195 | #Arrays |
---|
196 | #w_C = W.centroid_values |
---|
197 | #z_C = Z.centroid_values |
---|
198 | |
---|
199 | #num_min = num.min(w_C-z_C) |
---|
200 | #if num_min < 0.0: |
---|
201 | # print 'num.min(w_C-z_C)', num_min |
---|
202 | |
---|
203 | #print 'using swb2 distribute_to_vertices_and_edges' |
---|
204 | if self.use_edge_limiter: |
---|
205 | #distribute_using_edge_limiter(self) |
---|
206 | self.distribute_using_edge_limiter() |
---|
207 | else: |
---|
208 | #distribute_using_vertex_limiter(self) |
---|
209 | self.distribute_using_vertex_limiter() |
---|
210 | |
---|
211 | #def distribute_to_vertices_and_edges(self): |
---|
212 | # """Distribution from centroids to edges specific to the SWW eqn. |
---|
213 | |
---|
214 | # In addition, all conserved quantities get distributed as per either a |
---|
215 | # constant (order==1) or a piecewise linear function (order==2). |
---|
216 | # |
---|
217 | |
---|
218 | # Precondition: |
---|
219 | # All conserved quantities defined at centroids and bed elevation defined at |
---|
220 | # edges. |
---|
221 | # |
---|
222 | # Postcondition |
---|
223 | # Evolved quantities defined at vertices and edges |
---|
224 | # """ |
---|
225 | # from swb2_domain import protect_against_infinitesimal_and_negative_heights |
---|
226 | |
---|
227 | # if self._order_ == 1: |
---|
228 | # for name in [ 'stage', 'xmomentum', 'ymomentum' ]: |
---|
229 | # Q = self.quantities[name] |
---|
230 | # Q.extrapolate_first_order() |
---|
231 | # elif self._order_ == 2: |
---|
232 | # for name in [ 'stage', 'xmomentum', 'ymomentum' ]: |
---|
233 | # Q = self.quantities[name] |
---|
234 | # #Q.extrapolate_second_order_and_limit_by_edge() |
---|
235 | # Q.extrapolate_second_order_and_limit_by_vertex() |
---|
236 | |
---|
237 | # #self.extrapolate_second_order_sw() |
---|
238 | # else: |
---|
239 | # raise Exception('Unknown order: %s' % str(self._order_)) |
---|
240 | |
---|
241 | |
---|
242 | # # Update other quantities |
---|
243 | # protect_against_infinitesimal_and_negative_heights(self) |
---|
244 | |
---|
245 | # |
---|
246 | # # Compute edge values by interpolation |
---|
247 | # for name in ['stage', 'xmomentum', 'ymomentum']: |
---|
248 | # Q = self.quantities[name] |
---|
249 | # Q.interpolate_from_vertices_to_edges() |
---|
250 | |
---|
251 | |
---|
252 | def distribute_using_vertex_limiter(domain): |
---|
253 | """Distribution from centroids to vertices specific to the SWW equation. |
---|
254 | |
---|
255 | Precondition: |
---|
256 | All quantities defined at centroids and bed elevation defined at |
---|
257 | vertices. |
---|
258 | |
---|
259 | Postcondition |
---|
260 | Conserved quantities defined at vertices |
---|
261 | """ |
---|
262 | #from swb2_domain import protect_against_infinitesimal_and_negative_heights |
---|
263 | |
---|
264 | #print 'using swb2 distribute_using_vertex_limiter' |
---|
265 | # Remove very thin layers of water |
---|
266 | domain.protect_against_infinitesimal_and_negative_heights() |
---|
267 | |
---|
268 | # Extrapolate all conserved quantities |
---|
269 | if domain.optimised_gradient_limiter: |
---|
270 | # MH090605 if second order, |
---|
271 | # perform the extrapolation and limiting on |
---|
272 | # all of the conserved quantities |
---|
273 | |
---|
274 | if (domain._order_ == 1): |
---|
275 | for name in domain.conserved_quantities: |
---|
276 | Q = domain.quantities[name] |
---|
277 | Q.extrapolate_first_order() |
---|
278 | elif domain._order_ == 2: |
---|
279 | domain.extrapolate_second_order_sw() |
---|
280 | else: |
---|
281 | raise Exception('Unknown order') |
---|
282 | else: |
---|
283 | # Old code: |
---|
284 | for name in domain.conserved_quantities: |
---|
285 | Q = domain.quantities[name] |
---|
286 | |
---|
287 | if domain._order_ == 1: |
---|
288 | Q.extrapolate_first_order() |
---|
289 | elif domain._order_ == 2: |
---|
290 | Q.extrapolate_second_order_and_limit_by_vertex() |
---|
291 | else: |
---|
292 | raise Exception('Unknown order') |
---|
293 | |
---|
294 | # Take bed elevation into account when water heights are small |
---|
295 | #balance_deep_and_shallow(domain) |
---|
296 | |
---|
297 | # Compute edge values by interpolation |
---|
298 | for name in domain.conserved_quantities: |
---|
299 | Q = domain.quantities[name] |
---|
300 | Q.interpolate_from_vertices_to_edges() |
---|
301 | |
---|
302 | |
---|
303 | def distribute_using_edge_limiter(domain): |
---|
304 | """Distribution from centroids to edges specific to the SWW eqn. |
---|
305 | |
---|
306 | Precondition: |
---|
307 | All quantities defined at centroids and bed elevation defined at |
---|
308 | vertices. |
---|
309 | |
---|
310 | Postcondition |
---|
311 | Conserved quantities defined at vertices |
---|
312 | """ |
---|
313 | #from swb2_domain import protect_against_infinitesimal_and_negative_heights |
---|
314 | |
---|
315 | # Remove very thin layers of water |
---|
316 | domain.protect_against_infinitesimal_and_negative_heights() |
---|
317 | |
---|
318 | #print 'using swb2 distribute_using_vertex_limiter' |
---|
319 | for name in domain.conserved_quantities: |
---|
320 | Q = domain.quantities[name] |
---|
321 | if domain._order_ == 1: |
---|
322 | Q.extrapolate_first_order() |
---|
323 | elif domain._order_ == 2: |
---|
324 | Q.extrapolate_second_order_and_limit_by_edge() |
---|
325 | else: |
---|
326 | raise Exception('Unknown order') |
---|
327 | |
---|
328 | #balance_deep_and_shallow(domain) |
---|
329 | |
---|
330 | # Compute edge values by interpolation |
---|
331 | for name in domain.conserved_quantities: |
---|
332 | Q = domain.quantities[name] |
---|
333 | Q.interpolate_from_vertices_to_edges() |
---|
334 | |
---|
335 | |
---|
336 | def update_centroids_of_velocities_and_height(self): |
---|
337 | """Calculate the centroid values of velocities and height based |
---|
338 | on the values of the quantities stage and x and y momentum |
---|
339 | |
---|
340 | Assumes that stage and momentum are up to date |
---|
341 | |
---|
342 | Useful for kinematic viscosity calculations |
---|
343 | """ |
---|
344 | |
---|
345 | # For shallow water we need to update height xvelocity and yvelocity |
---|
346 | |
---|
347 | #Shortcuts |
---|
348 | W = self.quantities['stage'] |
---|
349 | UH = self.quantities['xmomentum'] |
---|
350 | VH = self.quantities['ymomentum'] |
---|
351 | H = self.quantities['height'] |
---|
352 | Z = self.quantities['elevation'] |
---|
353 | U = self.quantities['xvelocity'] |
---|
354 | V = self.quantities['yvelocity'] |
---|
355 | |
---|
356 | #print num.min(W.centroid_values) |
---|
357 | |
---|
358 | # Make sure boundary values of conserved quantites |
---|
359 | # are consistent with value of functions at centroids |
---|
360 | #self.distribute_to_vertices_and_edges() |
---|
361 | Z.set_boundary_values_from_edges() |
---|
362 | |
---|
363 | #W.set_boundary_values_from_edges() |
---|
364 | #UH.set_boundary_values_from_edges() |
---|
365 | #VH.set_boundary_values_from_edges() |
---|
366 | |
---|
367 | # Update height values |
---|
368 | H.set_values(W.centroid_values-Z.centroid_values, location='centroids') |
---|
369 | H.set_boundary_values( num.where(W.boundary_values-Z.boundary_values>=0, |
---|
370 | W.boundary_values-Z.boundary_values, 0.0)) |
---|
371 | |
---|
372 | #assert num.min(H.centroid_values) >= 0 |
---|
373 | #assert num.min(H.boundary_values) >= 0 |
---|
374 | |
---|
375 | #Aliases |
---|
376 | uh_C = UH.centroid_values |
---|
377 | vh_C = VH.centroid_values |
---|
378 | h_C = H.centroid_values |
---|
379 | |
---|
380 | uh_B = UH.boundary_values |
---|
381 | vh_B = VH.boundary_values |
---|
382 | h_B = H.boundary_values |
---|
383 | |
---|
384 | H0 = 1.0e-8 |
---|
385 | |
---|
386 | U.set_values(uh_C/(h_C + H0/h_C), location='centroids') |
---|
387 | V.set_values(vh_C/(h_C + H0/h_C), location='centroids') |
---|
388 | |
---|
389 | U.set_boundary_values(uh_B/(h_B + H0/h_B)) |
---|
390 | V.set_boundary_values(vh_B/(h_B + H0/h_B)) |
---|