1 | """Class Domain - 2D triangular domains for finite-volume computations of |
---|
2 | conservation laws. |
---|
3 | |
---|
4 | This is the base class for various domain models, such as: the Advection |
---|
5 | implementation is a simple algorithm, mainly for testing purposes, and |
---|
6 | the standard Shallow Water Wave domain (simply known as Domain) is the |
---|
7 | standard for realistic simulation. |
---|
8 | |
---|
9 | |
---|
10 | Copyright 2004 |
---|
11 | Ole Nielsen, Stephen Roberts, Duncan Gray |
---|
12 | Geoscience Australia |
---|
13 | """ |
---|
14 | |
---|
15 | import types |
---|
16 | from time import time as walltime |
---|
17 | |
---|
18 | from anuga.abstract_2d_finite_volumes.neighbour_mesh import Mesh |
---|
19 | from pmesh2domain import pmesh_to_domain |
---|
20 | from region import Set_region as region_set_region |
---|
21 | from anuga.geometry.polygon import inside_polygon |
---|
22 | from anuga.abstract_2d_finite_volumes.util import get_textual_float |
---|
23 | from quantity import Quantity |
---|
24 | import anuga.utilities.log as log |
---|
25 | |
---|
26 | import numpy as num |
---|
27 | |
---|
28 | class Generic_Domain: |
---|
29 | ''' |
---|
30 | Generic computational Domain constructor. |
---|
31 | # @param source Name of mesh file or coords of mesh vertices. |
---|
32 | # @param triangles Mesh connectivity (see mesh.py for more information). |
---|
33 | # @param boundary (see mesh.py for more information) |
---|
34 | # @param conserved_quantities List of names of quantities to be conserved. |
---|
35 | # @param other_quantities List of names of other quantities. |
---|
36 | # @param tagged_elements ?? |
---|
37 | # @param geo_reference ?? |
---|
38 | # @param use_inscribed_circle ?? |
---|
39 | # @param mesh_filename ?? |
---|
40 | # @param use_cache ?? |
---|
41 | # @param verbose True if this method is to be verbose. |
---|
42 | # @param full_send_dict ?? |
---|
43 | # @param ghost_recv_dict ?? |
---|
44 | # @param processor ?? |
---|
45 | # @param numproc ?? |
---|
46 | # @param number_of_full_nodes ?? |
---|
47 | # @param number_of_full_triangles ?? |
---|
48 | ''' |
---|
49 | |
---|
50 | def __init__(self, source=None, |
---|
51 | triangles=None, |
---|
52 | boundary=None, |
---|
53 | conserved_quantities=None, |
---|
54 | evolved_quantities=None, |
---|
55 | other_quantities=None, |
---|
56 | tagged_elements=None, |
---|
57 | geo_reference=None, |
---|
58 | use_inscribed_circle=False, |
---|
59 | mesh_filename=None, |
---|
60 | use_cache=False, |
---|
61 | verbose=False, |
---|
62 | full_send_dict=None, |
---|
63 | ghost_recv_dict=None, |
---|
64 | processor=0, |
---|
65 | numproc=1, |
---|
66 | number_of_full_nodes=None, |
---|
67 | number_of_full_triangles=None): |
---|
68 | |
---|
69 | """Instantiate generic computational Domain. |
---|
70 | |
---|
71 | Input: |
---|
72 | source: Either a mesh filename or coordinates of mesh vertices. |
---|
73 | If it is a filename values specified for triangles will |
---|
74 | be overridden. |
---|
75 | triangles: Mesh connectivity (see mesh.py for more information) |
---|
76 | boundary: See mesh.py for more information |
---|
77 | |
---|
78 | conserved_quantities: List of quantity names entering the |
---|
79 | conservation equations |
---|
80 | evolved_quantities: List of all quantities that evolve |
---|
81 | other_quantities: List of other quantity names |
---|
82 | |
---|
83 | tagged_elements: |
---|
84 | ... |
---|
85 | """ |
---|
86 | |
---|
87 | number_of_full_nodes=None |
---|
88 | number_of_full_triangles=None |
---|
89 | |
---|
90 | # Determine whether source is a mesh filename or coordinates |
---|
91 | if type(source) == types.StringType: |
---|
92 | mesh_filename = source |
---|
93 | else: |
---|
94 | coordinates = source |
---|
95 | |
---|
96 | # In case a filename has been specified, extract content |
---|
97 | if mesh_filename is not None: |
---|
98 | coordinates, triangles, boundary, vertex_quantity_dict, \ |
---|
99 | tagged_elements, geo_reference = \ |
---|
100 | pmesh_to_domain(file_name=mesh_filename, |
---|
101 | use_cache=use_cache, |
---|
102 | verbose=verbose) |
---|
103 | |
---|
104 | # Initialise underlying mesh structure |
---|
105 | self.mesh = Mesh(coordinates, triangles, |
---|
106 | boundary=boundary, |
---|
107 | tagged_elements=tagged_elements, |
---|
108 | geo_reference=geo_reference, |
---|
109 | use_inscribed_circle=use_inscribed_circle, |
---|
110 | number_of_full_nodes=number_of_full_nodes, |
---|
111 | number_of_full_triangles=number_of_full_triangles, |
---|
112 | verbose=verbose) |
---|
113 | |
---|
114 | # Expose Mesh attributes (FIXME: Maybe turn into methods) |
---|
115 | self.triangles = self.mesh.triangles |
---|
116 | self.centroid_coordinates = self.mesh.centroid_coordinates |
---|
117 | self.vertex_coordinates = self.mesh.vertex_coordinates |
---|
118 | self.boundary = self.mesh.boundary |
---|
119 | self.neighbours = self.mesh.neighbours |
---|
120 | self.surrogate_neighbours = self.mesh.surrogate_neighbours |
---|
121 | self.neighbour_edges = self.mesh.neighbour_edges |
---|
122 | self.normals = self.mesh.normals |
---|
123 | self.edgelengths = self.mesh.edgelengths |
---|
124 | self.radii = self.mesh.radii |
---|
125 | self.areas = self.mesh.areas |
---|
126 | |
---|
127 | self.number_of_boundaries = self.mesh.number_of_boundaries |
---|
128 | self.number_of_full_nodes = self.mesh.number_of_full_nodes |
---|
129 | self.number_of_full_triangles = self.mesh.number_of_full_triangles |
---|
130 | self.number_of_triangles_per_node = \ |
---|
131 | self.mesh.number_of_triangles_per_node |
---|
132 | |
---|
133 | self.vertex_value_indices = self.mesh.vertex_value_indices |
---|
134 | self.number_of_triangles = self.mesh.number_of_triangles |
---|
135 | |
---|
136 | self.geo_reference = self.mesh.geo_reference |
---|
137 | |
---|
138 | if verbose: log.critical('Initialising Domain') |
---|
139 | |
---|
140 | # List of quantity names entering the conservation equations |
---|
141 | if conserved_quantities is None: |
---|
142 | self.conserved_quantities = [] |
---|
143 | else: |
---|
144 | self.conserved_quantities = conserved_quantities |
---|
145 | |
---|
146 | if evolved_quantities is None: |
---|
147 | self.evolved_quantities = self.conserved_quantities |
---|
148 | else: |
---|
149 | self.evolved_quantities = evolved_quantities |
---|
150 | |
---|
151 | # List of other quantity names |
---|
152 | if other_quantities is None: |
---|
153 | self.other_quantities = [] |
---|
154 | else: |
---|
155 | self.other_quantities = other_quantities |
---|
156 | |
---|
157 | # Test that conserved_quantities are stored in the first entries of |
---|
158 | # evolved_quantities |
---|
159 | for i, quantity in enumerate(self.conserved_quantities): |
---|
160 | msg = 'The conserved quantities must be the first entries of ' |
---|
161 | msg += 'evolved_quantities' |
---|
162 | assert quantity == self.evolved_quantities[i], msg |
---|
163 | |
---|
164 | |
---|
165 | # Build dictionary of Quantity instances keyed by quantity names |
---|
166 | self.quantities = {} |
---|
167 | |
---|
168 | for name in self.evolved_quantities: |
---|
169 | self.quantities[name] = Quantity(self) |
---|
170 | for name in self.other_quantities: |
---|
171 | self.quantities[name] = Quantity(self) |
---|
172 | |
---|
173 | # Create an empty list for explicit forcing terms |
---|
174 | self.forcing_terms = [] |
---|
175 | |
---|
176 | # Setup the ghost cell communication |
---|
177 | if full_send_dict is None: |
---|
178 | self.full_send_dict = {} |
---|
179 | else: |
---|
180 | self.full_send_dict = full_send_dict |
---|
181 | |
---|
182 | # List of other quantity names |
---|
183 | if ghost_recv_dict is None: |
---|
184 | self.ghost_recv_dict = {} |
---|
185 | else: |
---|
186 | self.ghost_recv_dict = ghost_recv_dict |
---|
187 | |
---|
188 | self.processor = processor |
---|
189 | self.numproc = numproc |
---|
190 | |
---|
191 | # Setup Communication Buffers |
---|
192 | if verbose: log.critical('Domain: Set up communication buffers ' |
---|
193 | '(parallel)') |
---|
194 | self.nsys = len(self.conserved_quantities) |
---|
195 | for key in self.full_send_dict: |
---|
196 | buffer_shape = self.full_send_dict[key][0].shape[0] |
---|
197 | self.full_send_dict[key].append(num.zeros((buffer_shape, self.nsys), |
---|
198 | num.float)) |
---|
199 | |
---|
200 | for key in self.ghost_recv_dict: |
---|
201 | buffer_shape = self.ghost_recv_dict[key][0].shape[0] |
---|
202 | self.ghost_recv_dict[key].append( \ |
---|
203 | num.zeros((buffer_shape, self.nsys), |
---|
204 | num.float)) |
---|
205 | |
---|
206 | # Setup cell full flag |
---|
207 | # =1 for full |
---|
208 | # =0 for ghost |
---|
209 | N = len(self) #number_of_elements |
---|
210 | self.number_of_elements = N |
---|
211 | self.tri_full_flag = num.ones(N, num.int) |
---|
212 | for i in self.ghost_recv_dict.keys(): |
---|
213 | for id in self.ghost_recv_dict[i][0]: |
---|
214 | self.tri_full_flag[id] = 0 |
---|
215 | |
---|
216 | # Test the assumption that all full triangles are store before |
---|
217 | # the ghost triangles. |
---|
218 | if not num.allclose(self.tri_full_flag[:self.number_of_full_nodes], 1): |
---|
219 | if self.numproc>1: |
---|
220 | log.critical('WARNING: Not all full triangles are store before ' |
---|
221 | 'ghost triangles') |
---|
222 | |
---|
223 | # Defaults |
---|
224 | from anuga.config import max_smallsteps, beta_w, epsilon |
---|
225 | from anuga.config import CFL |
---|
226 | from anuga.config import timestepping_method |
---|
227 | from anuga.config import protect_against_isolated_degenerate_timesteps |
---|
228 | from anuga.config import default_order |
---|
229 | from anuga.config import max_timestep, min_timestep |
---|
230 | |
---|
231 | self.beta_w = beta_w |
---|
232 | self.epsilon = epsilon |
---|
233 | self.protect_against_isolated_degenerate_timesteps = \ |
---|
234 | protect_against_isolated_degenerate_timesteps |
---|
235 | |
---|
236 | |
---|
237 | self.centroid_transmissive_bc = False |
---|
238 | self.set_default_order(default_order) |
---|
239 | |
---|
240 | self.smallsteps = 0 |
---|
241 | self.max_smallsteps = max_smallsteps |
---|
242 | self.number_of_steps = 0 |
---|
243 | self.number_of_first_order_steps = 0 |
---|
244 | self.CFL = CFL |
---|
245 | self.set_timestepping_method(timestepping_method) |
---|
246 | self.set_beta(beta_w) |
---|
247 | self.set_evolve_max_timestep(max_timestep) |
---|
248 | self.set_evolve_min_timestep(min_timestep) |
---|
249 | self.boundary_map = None # Will be populated by set_boundary |
---|
250 | |
---|
251 | # Model time |
---|
252 | self.time = 0.0 |
---|
253 | self.finaltime = None |
---|
254 | self.recorded_min_timestep = self.recorded_max_timestep = 0.0 |
---|
255 | self.starttime = 0 # Physical starttime if any |
---|
256 | # (0 is 1 Jan 1970 00:00:00) |
---|
257 | self.timestep = 0.0 |
---|
258 | self.flux_timestep = 0.0 |
---|
259 | |
---|
260 | self.last_walltime = walltime() |
---|
261 | |
---|
262 | # Monitoring |
---|
263 | self.quantities_to_be_monitored = None |
---|
264 | self.monitor_polygon = None |
---|
265 | self.monitor_time_interval = None |
---|
266 | self.monitor_indices = None |
---|
267 | |
---|
268 | # Checkpointing and storage |
---|
269 | from anuga.config import default_datadir |
---|
270 | |
---|
271 | self.datadir = default_datadir |
---|
272 | self.simulation_name = 'domain' |
---|
273 | self.checkpoint = False |
---|
274 | |
---|
275 | # To avoid calculating the flux across each edge twice, keep an integer |
---|
276 | # (boolean) array, to be used during the flux calculation. |
---|
277 | N = len(self) # Number_of_triangles |
---|
278 | self.already_computed_flux = num.zeros((N, 3), num.int) |
---|
279 | |
---|
280 | # Storage for maximal speeds computed for each triangle by |
---|
281 | # compute_fluxes. |
---|
282 | # This is used for diagnostics only (reset at every yieldstep) |
---|
283 | self.max_speed = num.zeros(N, num.float) |
---|
284 | |
---|
285 | if mesh_filename is not None: |
---|
286 | # If the mesh file passed any quantity values, |
---|
287 | # initialise with these values. |
---|
288 | if verbose: log.critical('Domain: Initialising quantity values') |
---|
289 | self.set_quantity_vertices_dict(vertex_quantity_dict) |
---|
290 | |
---|
291 | if verbose: log.critical('Domain: Done') |
---|
292 | |
---|
293 | ###### |
---|
294 | # Expose underlying Mesh functionality |
---|
295 | ###### |
---|
296 | |
---|
297 | def __len__(self): |
---|
298 | return len(self.mesh) |
---|
299 | |
---|
300 | def get_centroid_coordinates(self, *args, **kwargs): |
---|
301 | return self.mesh.get_centroid_coordinates(*args, **kwargs) |
---|
302 | |
---|
303 | def get_radii(self, *args, **kwargs): |
---|
304 | return self.mesh.get_radii(*args, **kwargs) |
---|
305 | |
---|
306 | def get_areas(self, *args, **kwargs): |
---|
307 | return self.mesh.get_areas(*args, **kwargs) |
---|
308 | |
---|
309 | def get_area(self, *args, **kwargs): |
---|
310 | return self.mesh.get_area(*args, **kwargs) |
---|
311 | |
---|
312 | def get_vertex_coordinates(self, *args, **kwargs): |
---|
313 | return self.mesh.get_vertex_coordinates(*args, **kwargs) |
---|
314 | |
---|
315 | def get_vertex_coordinate(self, *args, **kwargs): |
---|
316 | return self.mesh.get_vertex_coordinate(*args, **kwargs) |
---|
317 | |
---|
318 | def get_edge_midpoint_coordinates(self, *args, **kwargs): |
---|
319 | return self.mesh.get_edge_midpoint_coordinates(*args, **kwargs) |
---|
320 | |
---|
321 | def get_edge_midpoint_coordinate(self, *args, **kwargs): |
---|
322 | return self.mesh.get_edge_midpoint_coordinate(*args, **kwargs) |
---|
323 | |
---|
324 | def get_triangles(self, *args, **kwargs): |
---|
325 | return self.mesh.get_triangles(*args, **kwargs) |
---|
326 | |
---|
327 | def get_nodes(self, *args, **kwargs): |
---|
328 | return self.mesh.get_nodes(*args, **kwargs) |
---|
329 | |
---|
330 | def get_number_of_nodes(self, *args, **kwargs): |
---|
331 | return self.mesh.get_number_of_nodes(*args, **kwargs) |
---|
332 | |
---|
333 | def get_number_of_triangles(self, *args, **kwargs): |
---|
334 | return self.mesh.get_number_of_triangles(*args, **kwargs) |
---|
335 | |
---|
336 | def get_normal(self, *args, **kwargs): |
---|
337 | return self.mesh.get_normal(*args, **kwargs) |
---|
338 | |
---|
339 | def get_triangle_containing_point(self, *args, **kwargs): |
---|
340 | return self.mesh.get_triangle_containing_point(*args, **kwargs) |
---|
341 | |
---|
342 | def get_intersecting_segments(self, *args, **kwargs): |
---|
343 | return self.mesh.get_intersecting_segments(*args, **kwargs) |
---|
344 | |
---|
345 | def get_disconnected_triangles(self, *args, **kwargs): |
---|
346 | return self.mesh.get_disconnected_triangles(*args, **kwargs) |
---|
347 | |
---|
348 | def get_boundary_tags(self, *args, **kwargs): |
---|
349 | return self.mesh.get_boundary_tags(*args, **kwargs) |
---|
350 | |
---|
351 | def get_boundary_polygon(self, *args, **kwargs): |
---|
352 | return self.mesh.get_boundary_polygon(*args, **kwargs) |
---|
353 | |
---|
354 | # FIXME(Ole): This doesn't seem to be required |
---|
355 | def get_number_of_triangles_per_node(self, *args, **kwargs): |
---|
356 | return self.mesh.get_number_of_triangles_per_node(*args, **kwargs) |
---|
357 | |
---|
358 | def get_triangles_and_vertices_per_node(self, *args, **kwargs): |
---|
359 | return self.mesh.get_triangles_and_vertices_per_node(*args, **kwargs) |
---|
360 | |
---|
361 | def get_interpolation_object(self, *args, **kwargs): |
---|
362 | return self.mesh.get_interpolation_object(*args, **kwargs) |
---|
363 | |
---|
364 | def get_tagged_elements(self, *args, **kwargs): |
---|
365 | return self.mesh.get_tagged_elements(*args, **kwargs) |
---|
366 | |
---|
367 | def get_lone_vertices(self, *args, **kwargs): |
---|
368 | return self.mesh.get_lone_vertices(*args, **kwargs) |
---|
369 | |
---|
370 | def get_unique_vertices(self, *args, **kwargs): |
---|
371 | return self.mesh.get_unique_vertices(*args, **kwargs) |
---|
372 | |
---|
373 | def get_georeference(self, *args, **kwargs): |
---|
374 | return self.mesh.get_georeference(*args, **kwargs) |
---|
375 | |
---|
376 | def set_georeference(self, *args, **kwargs): |
---|
377 | self.mesh.set_georeference(*args, **kwargs) |
---|
378 | |
---|
379 | def build_tagged_elements_dictionary(self, *args, **kwargs): |
---|
380 | self.mesh.build_tagged_elements_dictionary(*args, **kwargs) |
---|
381 | |
---|
382 | def statistics(self, *args, **kwargs): |
---|
383 | return self.mesh.statistics(*args, **kwargs) |
---|
384 | |
---|
385 | def get_extent(self, *args, **kwargs): |
---|
386 | return self.mesh.get_extent(*args, **kwargs) |
---|
387 | |
---|
388 | ## |
---|
389 | # @brief Get conserved quantities for a volume. |
---|
390 | # @param vol_id ID of the volume we want the conserved quantities for. |
---|
391 | # @param vertex If specified, use as index for edge values. |
---|
392 | # @param edge If specified, use as index for edge values. |
---|
393 | # @return Vector of conserved quantities. |
---|
394 | # @note If neither 'vertex' or 'edge' specified, use centroid values. |
---|
395 | # @note If both 'vertex' and 'edge' specified, raise exception. |
---|
396 | def get_conserved_quantities(self, vol_id, |
---|
397 | vertex=None, |
---|
398 | edge=None): |
---|
399 | """Get conserved quantities at volume vol_id. |
---|
400 | |
---|
401 | If vertex is specified use it as index for vertex values |
---|
402 | If edge is specified use it as index for edge values |
---|
403 | If neither are specified use centroid values |
---|
404 | If both are specified an exeception is raised |
---|
405 | |
---|
406 | Return value: Vector of length == number_of_conserved quantities |
---|
407 | """ |
---|
408 | |
---|
409 | if not (vertex is None or edge is None): |
---|
410 | msg = 'Values for both vertex and edge was specified.' |
---|
411 | msg += 'Only one (or none) is allowed.' |
---|
412 | raise Exception, msg |
---|
413 | |
---|
414 | q = num.zeros(len(self.conserved_quantities), num.float) |
---|
415 | |
---|
416 | for i, name in enumerate(self.conserved_quantities): |
---|
417 | Q = self.quantities[name] |
---|
418 | if vertex is not None: |
---|
419 | q[i] = Q.vertex_values[vol_id, vertex] |
---|
420 | elif edge is not None: |
---|
421 | q[i] = Q.edge_values[vol_id, edge] |
---|
422 | else: |
---|
423 | q[i] = Q.centroid_values[vol_id] |
---|
424 | |
---|
425 | return q |
---|
426 | |
---|
427 | ## |
---|
428 | # @brief Get evolved quantities for a volume. |
---|
429 | # @param vol_id ID of the volume we want the conserved quantities for. |
---|
430 | # @param vertex If specified, use as index for edge values. |
---|
431 | # @param edge If specified, use as index for edge values. |
---|
432 | # @return Vector of conserved quantities. |
---|
433 | # @note If neither 'vertex' or 'edge' specified, use centroid values. |
---|
434 | # @note If both 'vertex' and 'edge' specified, raise exception. |
---|
435 | def get_evolved_quantities(self, vol_id, |
---|
436 | vertex=None, |
---|
437 | edge=None): |
---|
438 | """Get evolved quantities at volume vol_id. |
---|
439 | |
---|
440 | If vertex is specified use it as index for vertex values |
---|
441 | If edge is specified use it as index for edge values |
---|
442 | If neither are specified use centroid values |
---|
443 | If both are specified an exeception is raised |
---|
444 | |
---|
445 | Return value: Vector of length == number_of_conserved quantities |
---|
446 | """ |
---|
447 | |
---|
448 | if not (vertex is None or edge is None): |
---|
449 | msg = 'Values for both vertex and edge was specified.' |
---|
450 | msg += 'Only one (or none) is allowed.' |
---|
451 | raise Exception, msg |
---|
452 | |
---|
453 | q = num.zeros(len(self.evolved_quantities), num.float) |
---|
454 | |
---|
455 | for i, name in enumerate(self.evolved_quantities): |
---|
456 | Q = self.quantities[name] |
---|
457 | if vertex is not None: |
---|
458 | q[i] = Q.vertex_values[vol_id, vertex] |
---|
459 | elif edge is not None: |
---|
460 | q[i] = Q.edge_values[vol_id, edge] |
---|
461 | else: |
---|
462 | q[i] = Q.centroid_values[vol_id] |
---|
463 | |
---|
464 | return q |
---|
465 | |
---|
466 | ## |
---|
467 | # @brief |
---|
468 | # @param flag |
---|
469 | def set_CFL(self, cfl=1.0): |
---|
470 | """Set CFL parameter, warn if greater than 1.0 |
---|
471 | """ |
---|
472 | if cfl > 1.0: |
---|
473 | self.CFL = cfl |
---|
474 | log.warn('Setting CFL > 1.0') |
---|
475 | |
---|
476 | assert cfl > 0.0 |
---|
477 | self.CFL = cfl |
---|
478 | |
---|
479 | |
---|
480 | |
---|
481 | ## |
---|
482 | # @brief Set the relative model time. |
---|
483 | # @param time The new model time (seconds). |
---|
484 | def set_time(self, time=0.0): |
---|
485 | """Set the model time (seconds).""" |
---|
486 | |
---|
487 | # FIXME: this is setting the relative time |
---|
488 | # Note that get_time and set_time are now not symmetric |
---|
489 | |
---|
490 | self.time = time |
---|
491 | |
---|
492 | ## |
---|
493 | # @brief Get the model time. |
---|
494 | # @return The absolute model time (seconds). |
---|
495 | def get_time(self): |
---|
496 | """Get the absolute model time (seconds).""" |
---|
497 | |
---|
498 | return self.time + self.starttime |
---|
499 | |
---|
500 | ## |
---|
501 | # @brief Set the default beta for limiting. |
---|
502 | # @param beta The new beta value. |
---|
503 | def set_beta(self, beta): |
---|
504 | """Set default beta for limiting.""" |
---|
505 | |
---|
506 | self.beta = beta |
---|
507 | for name in self.quantities: |
---|
508 | Q = self.quantities[name] |
---|
509 | Q.set_beta(beta) |
---|
510 | |
---|
511 | ## |
---|
512 | # @brief Get the beta value used for limiting. |
---|
513 | # @return The beta value used for limiting. |
---|
514 | def get_beta(self): |
---|
515 | """Get default beta for limiting.""" |
---|
516 | |
---|
517 | return self.beta |
---|
518 | |
---|
519 | |
---|
520 | ## |
---|
521 | # @brief Set the behaviour of the transmissive boundary condition |
---|
522 | # @param flag. True or False flag |
---|
523 | def set_centroid_transmissive_bc(self, flag): |
---|
524 | """Set behaviour of the transmissive boundary condition, namely |
---|
525 | calculate the BC using the centroid value of neighbouring cell |
---|
526 | or the calculated edge value. |
---|
527 | |
---|
528 | Centroid value is safer. |
---|
529 | |
---|
530 | Some of the limiters (extrapolate_second_order_and_limit_by_edge) |
---|
531 | don't limit boundary edge values (so that linear functions are reconstructed), |
---|
532 | |
---|
533 | In this case it is possible for a run away inflow to occur at a transmissive |
---|
534 | boundary. In this case set centroid_transmissive_bc to True""" |
---|
535 | |
---|
536 | self.centroid_transmissive_bc = flag |
---|
537 | |
---|
538 | ## |
---|
539 | # @brief Get the centroid_transmissive_bc flag |
---|
540 | # @return The beta value used for limiting. |
---|
541 | def get_centroid_transmissive_bc(self): |
---|
542 | """Get value of centroid_transmissive_bc flag.""" |
---|
543 | |
---|
544 | return self.centroid_transmissive_bc |
---|
545 | |
---|
546 | |
---|
547 | ## |
---|
548 | # @brief Set the max timestep for time evolution |
---|
549 | # @param max_timestep The new max timestep value. |
---|
550 | def set_evolve_max_timestep(self, max_timestep): |
---|
551 | """Set default max_timestep for evolving.""" |
---|
552 | |
---|
553 | self.evolve_max_timestep = max_timestep |
---|
554 | |
---|
555 | |
---|
556 | ## |
---|
557 | # @brief Get the max timestep for time evolution |
---|
558 | # @return The max timestep value. |
---|
559 | def get_evolve_max_timestep(self): |
---|
560 | """Set default max_timestep for evolving.""" |
---|
561 | |
---|
562 | return self.evolve_max_timestep |
---|
563 | |
---|
564 | ## |
---|
565 | # @brief Set the min timestep for time evolution |
---|
566 | # @param min_timestep The new min timestep value. |
---|
567 | def set_evolve_min_timestep(self, min_timestep): |
---|
568 | """Set default min_timestep for evolving.""" |
---|
569 | |
---|
570 | self.evolve_min_timestep = min_timestep |
---|
571 | |
---|
572 | |
---|
573 | ## |
---|
574 | # @brief Get the min timestep for time evolution |
---|
575 | # @return The min timestep value. |
---|
576 | def get_evolve_min_timestep(self): |
---|
577 | """Set default max_timestep for evolving.""" |
---|
578 | |
---|
579 | return self.evolve_min_timestep |
---|
580 | |
---|
581 | |
---|
582 | |
---|
583 | ## |
---|
584 | # @brief Set default (spatial) order. |
---|
585 | # @param n The new spatial order value. |
---|
586 | # @note If 'n' is not 1 or 2, raise exception. |
---|
587 | def set_default_order(self, n): |
---|
588 | """Set default (spatial) order to either 1 or 2.""" |
---|
589 | |
---|
590 | msg = 'Default order must be either 1 or 2. I got %s' % n |
---|
591 | assert n in [1,2], msg |
---|
592 | |
---|
593 | self.default_order = n |
---|
594 | self._order_ = self.default_order |
---|
595 | |
---|
596 | ## |
---|
597 | # @brief Set values of named quantities. |
---|
598 | # @param quantity_dict Dictionary containing name/value pairs. |
---|
599 | def set_quantity_vertices_dict(self, quantity_dict): |
---|
600 | """Set values for named quantities. |
---|
601 | Supplied dictionary contains name/value pairs: |
---|
602 | |
---|
603 | name: Name of quantity |
---|
604 | value: Compatible list, numeric array, const or function (see below) |
---|
605 | |
---|
606 | The values will be stored in elements following their internal ordering. |
---|
607 | """ |
---|
608 | |
---|
609 | # FIXME: Could we name this a bit more intuitively |
---|
610 | # E.g. set_quantities_from_dictionary |
---|
611 | for key in quantity_dict.keys(): |
---|
612 | self.set_quantity(key, quantity_dict[key], location='vertices') |
---|
613 | |
---|
614 | ## |
---|
615 | # @brief Set value(s) for a named quantity. |
---|
616 | # @param name Name of quantity to be updated. |
---|
617 | # @param args Positional args. |
---|
618 | # @param kwargs Keyword args. |
---|
619 | # @note If 'kwargs' dict has 'expression' key, evaluate expression. |
---|
620 | def set_quantity(self, name, |
---|
621 | *args, **kwargs): |
---|
622 | """Set values for named quantity |
---|
623 | |
---|
624 | One keyword argument is documented here: |
---|
625 | expression = None, # Arbitrary expression |
---|
626 | |
---|
627 | expression: |
---|
628 | Arbitrary expression involving quantity names |
---|
629 | |
---|
630 | See Quantity.set_values for further documentation. |
---|
631 | """ |
---|
632 | |
---|
633 | # Do the expression stuff |
---|
634 | if kwargs.has_key('expression'): |
---|
635 | expression = kwargs['expression'] |
---|
636 | del kwargs['expression'] |
---|
637 | |
---|
638 | Q = self.create_quantity_from_expression(expression) |
---|
639 | kwargs['quantity'] = Q |
---|
640 | |
---|
641 | # Assign values |
---|
642 | self.quantities[name].set_values(*args, **kwargs) |
---|
643 | |
---|
644 | ## |
---|
645 | # @brief Add to a named quantity value. |
---|
646 | # @param name Name of quantity to be added to. |
---|
647 | # @param args Positional args. |
---|
648 | # @param kwargs Keyword args. |
---|
649 | # @note If 'kwargs' dict has 'expression' key, evaluate expression. |
---|
650 | def add_quantity(self, name, |
---|
651 | *args, **kwargs): |
---|
652 | """Add values to a named quantity |
---|
653 | |
---|
654 | E.g add_quantity('elevation', X) |
---|
655 | |
---|
656 | Option are the same as in set_quantity. |
---|
657 | """ |
---|
658 | |
---|
659 | # Do the expression stuff |
---|
660 | if kwargs.has_key('expression'): |
---|
661 | expression = kwargs['expression'] |
---|
662 | Q2 = self.create_quantity_from_expression(expression) |
---|
663 | else: |
---|
664 | # Create new temporary quantity |
---|
665 | Q2 = Quantity(self) |
---|
666 | |
---|
667 | # Assign specified values to temporary quantity |
---|
668 | Q2.set_values(*args, **kwargs) |
---|
669 | |
---|
670 | # Add temporary quantity to named quantity |
---|
671 | Q1 = self.get_quantity(name) |
---|
672 | self.set_quantity(name, Q1 + Q2) |
---|
673 | |
---|
674 | ## |
---|
675 | # @brief Get list of quantity names for the Domain. |
---|
676 | # @return List of quantity names. |
---|
677 | def get_quantity_names(self): |
---|
678 | """Get a list of all the quantity names that this domain is aware of. |
---|
679 | Any value in the result should be a valid input to get_quantity. |
---|
680 | """ |
---|
681 | |
---|
682 | return self.quantities.keys() |
---|
683 | |
---|
684 | ## |
---|
685 | # @brief Get a quantity object. |
---|
686 | # @param name Name of the quantity value. |
---|
687 | # @param location ?? |
---|
688 | # @param indices ?? |
---|
689 | # @return The quantity value object. |
---|
690 | # @note 'location' and 'indices' are unused. |
---|
691 | def get_quantity(self, name, |
---|
692 | location='vertices', |
---|
693 | indices = None): |
---|
694 | """Get pointer to quantity object. |
---|
695 | |
---|
696 | name: Name of quantity |
---|
697 | |
---|
698 | See methods inside the quantity object for more options |
---|
699 | |
---|
700 | FIXME: clean input args |
---|
701 | """ |
---|
702 | |
---|
703 | return self.quantities[name] #.get_values( location, indices = indices) |
---|
704 | |
---|
705 | ## |
---|
706 | # @brief Create a quantity value from an expression. |
---|
707 | # @param expression The expression (string) to be evaluated. |
---|
708 | # @return The expression value, evaluated from this Domain's quantities. |
---|
709 | # @note Valid expression operators are as defined in class Quantity. |
---|
710 | def create_quantity_from_expression(self, expression): |
---|
711 | """Create new quantity from other quantities using arbitrary expression. |
---|
712 | |
---|
713 | Combine existing quantities in domain using expression and return |
---|
714 | result as a new quantity. |
---|
715 | |
---|
716 | Note, the new quantity could e.g. be used in set_quantity |
---|
717 | |
---|
718 | Valid expressions are limited to operators defined in class Quantity |
---|
719 | |
---|
720 | Examples creating derived quantities: |
---|
721 | Depth = domain.create_quantity_from_expression('stage-elevation') |
---|
722 | exp = '(xmomentum*xmomentum + ymomentum*ymomentum)**0.5' |
---|
723 | Absolute_momentum = domain.create_quantity_from_expression(exp) |
---|
724 | """ |
---|
725 | |
---|
726 | from anuga.abstract_2d_finite_volumes.util import\ |
---|
727 | apply_expression_to_dictionary |
---|
728 | |
---|
729 | return apply_expression_to_dictionary(expression, self.quantities) |
---|
730 | |
---|
731 | ## |
---|
732 | # @brief Associate boundary objects with tagged boundary segments. |
---|
733 | # @param boundary_map A dict of boundary objects keyed by symbolic tags to |
---|
734 | # matched against tags in the internal dictionary |
---|
735 | # self.boundary. |
---|
736 | def set_boundary(self, boundary_map): |
---|
737 | """Associate boundary objects with tagged boundary segments. |
---|
738 | |
---|
739 | Input boundary_map is a dictionary of boundary objects keyed |
---|
740 | by symbolic tags to matched against tags in the internal dictionary |
---|
741 | self.boundary. |
---|
742 | |
---|
743 | As result one pointer to a boundary object is stored for each vertex |
---|
744 | in the list self.boundary_objects. |
---|
745 | More entries may point to the same boundary object |
---|
746 | |
---|
747 | Schematically the mapping is from two dictionaries to one list |
---|
748 | where the index is used as pointer to the boundary_values arrays |
---|
749 | within each quantity. |
---|
750 | |
---|
751 | self.boundary: (vol_id, edge_id): tag |
---|
752 | boundary_map (input): tag: boundary_object |
---|
753 | ---------------------------------------------- |
---|
754 | self.boundary_objects: ((vol_id, edge_id), boundary_object) |
---|
755 | |
---|
756 | Pre-condition: |
---|
757 | self.boundary has been built. |
---|
758 | |
---|
759 | Post-condition: |
---|
760 | self.boundary_objects is built |
---|
761 | |
---|
762 | If a tag from the domain doesn't appear in the input dictionary an |
---|
763 | exception is raised. |
---|
764 | However, if a tag is not used to the domain, no error is thrown. |
---|
765 | FIXME: This would lead to implementation of a default boundary condition |
---|
766 | |
---|
767 | Note: If a segment is listed in the boundary dictionary and if it is |
---|
768 | not None, it *will* become a boundary - even if there is a neighbouring |
---|
769 | triangle. This would be the case for internal boundaries. |
---|
770 | |
---|
771 | Boundary objects that are None will be skipped. |
---|
772 | |
---|
773 | If a boundary_map has already been set (i.e. set_boundary has been |
---|
774 | called before), the old boundary map will be updated with new values. |
---|
775 | The new map need not define all boundary tags, and can thus change only |
---|
776 | those that are needed. |
---|
777 | |
---|
778 | FIXME: If set_boundary is called multiple times and if Boundary |
---|
779 | object is changed into None, the neighbour structure will not be |
---|
780 | restored!!! |
---|
781 | """ |
---|
782 | |
---|
783 | if self.boundary_map is None: |
---|
784 | # This the first call to set_boundary. Store |
---|
785 | # map for later updates and for use with boundary_stats. |
---|
786 | self.boundary_map = boundary_map |
---|
787 | else: |
---|
788 | # This is a modification of an already existing map |
---|
789 | # Update map an proceed normally |
---|
790 | for key in boundary_map.keys(): |
---|
791 | self.boundary_map[key] = boundary_map[key] |
---|
792 | |
---|
793 | # FIXME (Ole): Try to remove the sorting and fix test_mesh.py |
---|
794 | x = self.boundary.keys() |
---|
795 | x.sort() |
---|
796 | |
---|
797 | # Loop through edges that lie on the boundary and associate them with |
---|
798 | # callable boundary objects depending on their tags |
---|
799 | self.boundary_objects = [] |
---|
800 | for k, (vol_id, edge_id) in enumerate(x): |
---|
801 | tag = self.boundary[(vol_id, edge_id)] |
---|
802 | |
---|
803 | if self.boundary_map.has_key(tag): |
---|
804 | B = self.boundary_map[tag] # Get callable boundary object |
---|
805 | |
---|
806 | if B is not None: |
---|
807 | self.boundary_objects.append(((vol_id, edge_id), B)) |
---|
808 | self.neighbours[vol_id, edge_id] = \ |
---|
809 | -len(self.boundary_objects) |
---|
810 | else: |
---|
811 | pass |
---|
812 | #FIXME: Check and perhaps fix neighbour structure |
---|
813 | else: |
---|
814 | msg = 'ERROR (domain.py): Tag "%s" has not been ' %tag |
---|
815 | msg += 'bound to a boundary object.\n' |
---|
816 | msg += 'All boundary tags defined in domain must appear ' |
---|
817 | msg += 'in set_boundary.\n' |
---|
818 | msg += 'The tags are: %s' %self.get_boundary_tags() |
---|
819 | raise Exception, msg |
---|
820 | |
---|
821 | ## |
---|
822 | # @brief Set quantities based on a regional tag. |
---|
823 | # @param args |
---|
824 | # @param kwargs |
---|
825 | def set_region(self, *args, **kwargs): |
---|
826 | """Set quantities based on a regional tag. |
---|
827 | |
---|
828 | It is most often called with the following parameters; |
---|
829 | (self, tag, quantity, X, location='vertices') |
---|
830 | tag: the name of the regional tag used to specify the region |
---|
831 | quantity: Name of quantity to change |
---|
832 | X: const or function - how the quantity is changed |
---|
833 | location: Where values are to be stored. |
---|
834 | Permissible options are: vertices, centroid and unique vertices |
---|
835 | |
---|
836 | A callable region class or a list of callable region classes |
---|
837 | can also be passed into this function. |
---|
838 | """ |
---|
839 | |
---|
840 | if len(args) == 1: |
---|
841 | self._set_region(*args, **kwargs) |
---|
842 | else: |
---|
843 | # Assume it is arguments for the region.set_region function |
---|
844 | func = region_set_region(*args, **kwargs) |
---|
845 | self._set_region(func) |
---|
846 | |
---|
847 | ## |
---|
848 | # @brief ?? |
---|
849 | # @param functions A list or tuple of ?? |
---|
850 | def _set_region(self, functions): |
---|
851 | # coerce to an iterable (list or tuple) |
---|
852 | if type(functions) not in [types.ListType, types.TupleType]: |
---|
853 | functions = [functions] |
---|
854 | |
---|
855 | # The order of functions in the list is used. |
---|
856 | tagged_elements = self.get_tagged_elements() |
---|
857 | for function in functions: |
---|
858 | for tag in tagged_elements.keys(): |
---|
859 | function(tag, tagged_elements[tag], self) |
---|
860 | |
---|
861 | ## |
---|
862 | # @brief Specify the quantities which will be monitored for extrema. |
---|
863 | # @param q Single or list of quantity names to monitor. |
---|
864 | # @param polygon If specified, monitor only triangles inside polygon. |
---|
865 | # @param time_interval If specified, monitor only timesteps inside interval. |
---|
866 | # @note If 'q' is None, do no monitoring. |
---|
867 | def set_quantities_to_be_monitored(self, q, |
---|
868 | polygon=None, |
---|
869 | time_interval=None): |
---|
870 | """Specify which quantities will be monitored for extrema. |
---|
871 | |
---|
872 | q must be either: |
---|
873 | - the name of a quantity or derived quantity such as 'stage-elevation' |
---|
874 | - a list of quantity names |
---|
875 | - None |
---|
876 | |
---|
877 | In the two first cases, the named quantities will be monitored at |
---|
878 | each internal timestep |
---|
879 | |
---|
880 | If q is None, monitoring will be switched off altogether. |
---|
881 | |
---|
882 | polygon (if specified) will only monitor triangles inside polygon. |
---|
883 | If omitted all triangles will be included. |
---|
884 | |
---|
885 | time_interval, if specified, will restrict monitoring to time steps in |
---|
886 | that interval. If omitted all timesteps will be included. |
---|
887 | """ |
---|
888 | |
---|
889 | from anuga.abstract_2d_finite_volumes.util import\ |
---|
890 | apply_expression_to_dictionary |
---|
891 | |
---|
892 | if q is None: |
---|
893 | self.quantities_to_be_monitored = None |
---|
894 | self.monitor_polygon = None |
---|
895 | self.monitor_time_interval = None |
---|
896 | self.monitor_indices = None |
---|
897 | return |
---|
898 | |
---|
899 | # coerce 'q' to a list if it's a string |
---|
900 | if isinstance(q, basestring): |
---|
901 | q = [q] |
---|
902 | |
---|
903 | # Check correctness and initialise |
---|
904 | self.quantities_to_be_monitored = {} |
---|
905 | for quantity_name in q: |
---|
906 | msg = 'Quantity %s is not a valid conserved quantity' \ |
---|
907 | % quantity_name |
---|
908 | |
---|
909 | if not quantity_name in self.quantities: |
---|
910 | # See if this expression is valid |
---|
911 | apply_expression_to_dictionary(quantity_name, self.quantities) |
---|
912 | |
---|
913 | # Initialise extrema information |
---|
914 | info_block = {'min': None, # Min value |
---|
915 | 'max': None, # Max value |
---|
916 | 'min_location': None, # Argmin (x, y) |
---|
917 | 'max_location': None, # Argmax (x, y) |
---|
918 | 'min_time': None, # Argmin (t) |
---|
919 | 'max_time': None} # Argmax (t) |
---|
920 | |
---|
921 | self.quantities_to_be_monitored[quantity_name] = info_block |
---|
922 | |
---|
923 | if polygon is not None: |
---|
924 | # Check input |
---|
925 | if isinstance(polygon, basestring): |
---|
926 | # Check if multiple quantities were accidentally |
---|
927 | # given as separate argument rather than a list. |
---|
928 | msg = ('Multiple quantities must be specified in a list. ' |
---|
929 | 'Not as multiple arguments. ' |
---|
930 | 'I got "%s" as a second argument') % polygon |
---|
931 | |
---|
932 | if polygon in self.quantities: |
---|
933 | raise Exception, msg |
---|
934 | |
---|
935 | try: |
---|
936 | apply_expression_to_dictionary(polygon, self.quantities) |
---|
937 | except: |
---|
938 | # At least polygon wasn't expression involving quantitites |
---|
939 | pass |
---|
940 | else: |
---|
941 | raise Exception, msg |
---|
942 | |
---|
943 | # In any case, we don't allow polygon to be a string |
---|
944 | msg = ('argument "polygon" must not be a string: ' |
---|
945 | 'I got polygon="%s"') % polygon |
---|
946 | raise Exception, msg |
---|
947 | |
---|
948 | # Get indices for centroids that are inside polygon |
---|
949 | points = self.get_centroid_coordinates(absolute=True) |
---|
950 | self.monitor_indices = inside_polygon(points, polygon) |
---|
951 | |
---|
952 | if time_interval is not None: |
---|
953 | assert len(time_interval) == 2 |
---|
954 | |
---|
955 | self.monitor_polygon = polygon |
---|
956 | self.monitor_time_interval = time_interval |
---|
957 | |
---|
958 | ## |
---|
959 | # @brief Check Domain integrity. |
---|
960 | # @note Raises an exception if integrity breached. |
---|
961 | def check_integrity(self): |
---|
962 | self.mesh.check_integrity() |
---|
963 | |
---|
964 | for quantity in self.conserved_quantities: |
---|
965 | msg = 'Conserved quantities must be a subset of all quantities' |
---|
966 | assert quantity in self.quantities, msg |
---|
967 | |
---|
968 | |
---|
969 | for i, quantity in enumerate(self.conserved_quantities): |
---|
970 | msg = 'Conserved quantities must be the first entries ' |
---|
971 | msg += 'of evolved_quantities' |
---|
972 | assert quantity == self.evolved_quantities[i], msg |
---|
973 | |
---|
974 | |
---|
975 | ## |
---|
976 | # @brief Print timestep stats to stdout. |
---|
977 | # @param track_speeds If True, print smallest track speed. |
---|
978 | def write_time(self, track_speeds=False): |
---|
979 | log.critical(self.timestepping_statistics(track_speeds)) |
---|
980 | |
---|
981 | ## |
---|
982 | # @brief Get timestepping stats string. |
---|
983 | # @param track_speeds If True, report location of smallest timestep. |
---|
984 | # @param triangle_id If specified, use specific triangle. |
---|
985 | # @return A string containing timestep stats. |
---|
986 | def timestepping_statistics(self, track_speeds=False, |
---|
987 | triangle_id=None): |
---|
988 | """Return string with time stepping statistics |
---|
989 | |
---|
990 | Optional boolean keyword track_speeds decides whether to report |
---|
991 | location of smallest timestep as well as a histogram and percentile |
---|
992 | report. |
---|
993 | |
---|
994 | Optional keyword triangle_id can be used to specify a particular |
---|
995 | triangle rather than the one with the largest speed. |
---|
996 | """ |
---|
997 | |
---|
998 | from anuga.utilities.numerical_tools import histogram, create_bins |
---|
999 | |
---|
1000 | # qwidth determines the the width of the text field used for quantities |
---|
1001 | qwidth = self.qwidth = 12 |
---|
1002 | |
---|
1003 | msg = '' |
---|
1004 | |
---|
1005 | model_time = self.get_time() |
---|
1006 | if self.recorded_min_timestep == self.recorded_max_timestep: |
---|
1007 | msg += 'Time = %.4f, delta t = %.8f, steps=%d' \ |
---|
1008 | % (model_time, self.recorded_min_timestep, \ |
---|
1009 | self.number_of_steps) |
---|
1010 | elif self.recorded_min_timestep > self.recorded_max_timestep: |
---|
1011 | msg += 'Time = %.4f, steps=%d' \ |
---|
1012 | % (model_time, self.number_of_steps) |
---|
1013 | else: |
---|
1014 | msg += 'Time = %.4f, delta t in [%.8f, %.8f], steps=%d' \ |
---|
1015 | % (model_time, self.recorded_min_timestep, |
---|
1016 | self.recorded_max_timestep, self.number_of_steps) |
---|
1017 | |
---|
1018 | msg += ' (%ds)' % (walltime() - self.last_walltime) |
---|
1019 | self.last_walltime = walltime() |
---|
1020 | |
---|
1021 | if track_speeds is True: |
---|
1022 | msg += '\n' |
---|
1023 | |
---|
1024 | # Setup 10 bins for speed histogram |
---|
1025 | bins = create_bins(self.max_speed, 10) |
---|
1026 | hist = histogram(self.max_speed, bins) |
---|
1027 | |
---|
1028 | msg += '------------------------------------------------\n' |
---|
1029 | msg += ' Speeds in [%f, %f]\n' % (num.min(self.max_speed), |
---|
1030 | num.max(self.max_speed)) |
---|
1031 | msg += ' Histogram:\n' |
---|
1032 | |
---|
1033 | hi = bins[0] |
---|
1034 | for i, count in enumerate(hist): |
---|
1035 | lo = hi |
---|
1036 | if i+1 < len(bins): |
---|
1037 | # Open upper interval |
---|
1038 | hi = bins[i+1] |
---|
1039 | msg += ' [%f, %f[: %d\n' % (lo, hi, count) |
---|
1040 | else: |
---|
1041 | # Closed upper interval |
---|
1042 | hi = num.max(self.max_speed) |
---|
1043 | msg += ' [%f, %f]: %d\n' % (lo, hi, count) |
---|
1044 | |
---|
1045 | N = len(self.max_speed.flat) |
---|
1046 | if N > 10: |
---|
1047 | msg += ' Percentiles (10%):\n' |
---|
1048 | speed = self.max_speed.tolist() |
---|
1049 | speed.sort() |
---|
1050 | |
---|
1051 | k = 0 |
---|
1052 | lower = min(speed) |
---|
1053 | for i, a in enumerate(speed): |
---|
1054 | if i % (N/10) == 0 and i != 0: |
---|
1055 | # For every 10% of the sorted speeds |
---|
1056 | msg += ' %d speeds in [%f, %f]\n' % (i-k, lower, a) |
---|
1057 | lower = a |
---|
1058 | k = i |
---|
1059 | |
---|
1060 | msg += ' %d speeds in [%f, %f]\n'\ |
---|
1061 | % (N-k, lower, max(speed)) |
---|
1062 | |
---|
1063 | # Find index of largest computed flux speed |
---|
1064 | if triangle_id is None: |
---|
1065 | k = self.k = num.argmax(self.max_speed) |
---|
1066 | else: |
---|
1067 | errmsg = 'Triangle_id %d does not exist in mesh: %s' \ |
---|
1068 | % (triangle_id, str(self)) |
---|
1069 | assert 0 <= triangle_id < len(self), errmsg |
---|
1070 | k = self.k = triangle_id |
---|
1071 | |
---|
1072 | x, y = self.get_centroid_coordinates(absolute=True)[k] |
---|
1073 | radius = self.get_radii()[k] |
---|
1074 | area = self.get_areas()[k] |
---|
1075 | max_speed = self.max_speed[k] |
---|
1076 | |
---|
1077 | msg += ' Triangle #%d with centroid (%.4f, %.4f), ' % (k, x, y) |
---|
1078 | msg += 'area = %.4f and radius = %.4f ' % (area, radius) |
---|
1079 | if triangle_id is None: |
---|
1080 | msg += 'had the largest computed speed: %.6f m/s ' % (max_speed) |
---|
1081 | else: |
---|
1082 | msg += 'had computed speed: %.6f m/s ' % (max_speed) |
---|
1083 | |
---|
1084 | if max_speed > 0.0: |
---|
1085 | msg += '(timestep=%.6f)\n' % (radius/max_speed) |
---|
1086 | else: |
---|
1087 | msg += '(timestep=%.6f)\n' % (0) |
---|
1088 | |
---|
1089 | # Report all quantity values at vertices, edges and centroid |
---|
1090 | msg += ' Quantity' |
---|
1091 | msg += '------------\n' |
---|
1092 | for name in self.quantities: |
---|
1093 | q = self.quantities[name] |
---|
1094 | |
---|
1095 | V = q.get_values(location='vertices', indices=[k])[0] |
---|
1096 | E = q.get_values(location='edges', indices=[k])[0] |
---|
1097 | C = q.get_values(location='centroids', indices=[k]) |
---|
1098 | |
---|
1099 | s = ' %s: vertex_values = %.4f,\t %.4f,\t %.4f\n' \ |
---|
1100 | % (name.ljust(qwidth), V[0], V[1], V[2]) |
---|
1101 | |
---|
1102 | s += ' %s: edge_values = %.4f,\t %.4f,\t %.4f\n' \ |
---|
1103 | % (name.ljust(qwidth), E[0], E[1], E[2]) |
---|
1104 | |
---|
1105 | s += ' %s: centroid_value = %.4f\n' \ |
---|
1106 | % (name.ljust(qwidth), C[0]) |
---|
1107 | |
---|
1108 | msg += s |
---|
1109 | |
---|
1110 | return msg |
---|
1111 | |
---|
1112 | ## |
---|
1113 | # @brief Print boundary forcing stats at each timestep to stdout. |
---|
1114 | # @param quantities A name or list of names of quantities to report. |
---|
1115 | # @param tags A name or list of names of tags to report. |
---|
1116 | def write_boundary_statistics(self, quantities=None, tags=None): |
---|
1117 | log.critical(self.boundary_statistics(quantities, tags)) |
---|
1118 | |
---|
1119 | # @brief Get a string containing boundary forcing stats at each timestep. |
---|
1120 | # @param quantities A name or list of names of quantities to report. |
---|
1121 | # @param tags A name or list of names of tags to report. |
---|
1122 | # @note If 'quantities' is None, report all. Same for 'tags'. |
---|
1123 | def boundary_statistics(self, quantities=None, |
---|
1124 | tags=None): |
---|
1125 | """Output statistics about boundary forcing at each timestep |
---|
1126 | |
---|
1127 | Input: |
---|
1128 | quantities: either None, a string or a list of strings naming the |
---|
1129 | quantities to be reported |
---|
1130 | tags: either None, a string or a list of strings naming the |
---|
1131 | tags to be reported |
---|
1132 | |
---|
1133 | Example output: |
---|
1134 | Tag 'wall': |
---|
1135 | stage in [2, 5.5] |
---|
1136 | xmomentum in [] |
---|
1137 | ymomentum in [] |
---|
1138 | Tag 'ocean' |
---|
1139 | |
---|
1140 | If quantities are specified only report on those. Otherwise take all |
---|
1141 | conserved quantities. |
---|
1142 | If tags are specified only report on those, otherwise take all tags. |
---|
1143 | """ |
---|
1144 | |
---|
1145 | import types, string |
---|
1146 | |
---|
1147 | # Input checks |
---|
1148 | if quantities is None: |
---|
1149 | quantities = self.evolved_quantities |
---|
1150 | elif type(quantities) == types.StringType: |
---|
1151 | quantities = [quantities] #Turn it into a list |
---|
1152 | |
---|
1153 | msg = ('Keyword argument quantities must be either None, ' |
---|
1154 | 'string or list. I got %s') % str(quantities) |
---|
1155 | assert type(quantities) == types.ListType, msg |
---|
1156 | |
---|
1157 | if tags is None: |
---|
1158 | tags = self.get_boundary_tags() |
---|
1159 | elif type(tags) == types.StringType: |
---|
1160 | tags = [tags] #Turn it into a list |
---|
1161 | |
---|
1162 | msg = ('Keyword argument tags must be either None, ' |
---|
1163 | 'string or list. I got %s') % str(tags) |
---|
1164 | assert type(tags) == types.ListType, msg |
---|
1165 | |
---|
1166 | # Determine width of longest quantity name (for cosmetic purposes) |
---|
1167 | maxwidth = 0 |
---|
1168 | for name in quantities: |
---|
1169 | w = len(name) |
---|
1170 | if w > maxwidth: |
---|
1171 | maxwidth = w |
---|
1172 | |
---|
1173 | # Output statistics |
---|
1174 | msg = 'Boundary values at time %.4f:\n' % self.get_time() |
---|
1175 | for tag in tags: |
---|
1176 | msg += ' %s:\n' % tag |
---|
1177 | |
---|
1178 | for name in quantities: |
---|
1179 | q = self.quantities[name] |
---|
1180 | |
---|
1181 | # Find range of boundary values for tag and q |
---|
1182 | maxval = minval = None |
---|
1183 | for i, ((vol_id,edge_id),B) in enumerate(self.boundary_objects): |
---|
1184 | if self.boundary[(vol_id, edge_id)] == tag: |
---|
1185 | v = q.boundary_values[i] |
---|
1186 | if minval is None or v < minval: minval = v |
---|
1187 | if maxval is None or v > maxval: maxval = v |
---|
1188 | |
---|
1189 | if minval is None or maxval is None: |
---|
1190 | msg += (' Sorry no information available about' |
---|
1191 | ' tag %s and quantity %s\n') % (tag, name) |
---|
1192 | else: |
---|
1193 | msg += ' %s in [%12.8f, %12.8f]\n' \ |
---|
1194 | % (string.ljust(name, maxwidth), minval, maxval) |
---|
1195 | |
---|
1196 | return msg |
---|
1197 | |
---|
1198 | ## |
---|
1199 | # @brief Update extrema if requested by set_quantities_to_be_monitored. |
---|
1200 | def update_extrema(self): |
---|
1201 | """Update extrema if requested by set_quantities_to_be_monitored. |
---|
1202 | This data is used for reporting e.g. by running |
---|
1203 | print domain.quantity_statistics() |
---|
1204 | and may also stored in output files (see data_manager in shallow_water) |
---|
1205 | """ |
---|
1206 | |
---|
1207 | # Define a tolerance for extremum computations |
---|
1208 | from anuga.config import single_precision as epsilon |
---|
1209 | |
---|
1210 | if self.quantities_to_be_monitored is None: |
---|
1211 | return |
---|
1212 | |
---|
1213 | # Observe time interval restriction if any |
---|
1214 | if self.monitor_time_interval is not None and\ |
---|
1215 | (self.time < self.monitor_time_interval[0] or\ |
---|
1216 | self.time > self.monitor_time_interval[1]): |
---|
1217 | return |
---|
1218 | |
---|
1219 | # Update extrema for each specified quantity subject to |
---|
1220 | # polygon restriction (via monitor_indices). |
---|
1221 | for quantity_name in self.quantities_to_be_monitored: |
---|
1222 | |
---|
1223 | if quantity_name in self.quantities: |
---|
1224 | Q = self.get_quantity(quantity_name) |
---|
1225 | else: |
---|
1226 | Q = self.create_quantity_from_expression(quantity_name) |
---|
1227 | |
---|
1228 | info_block = self.quantities_to_be_monitored[quantity_name] |
---|
1229 | |
---|
1230 | # Update maximum |
---|
1231 | # (n > None is always True, but we check explicitly because |
---|
1232 | # of the epsilon) |
---|
1233 | maxval = Q.get_maximum_value(self.monitor_indices) |
---|
1234 | if info_block['max'] is None or \ |
---|
1235 | maxval > info_block['max'] + epsilon: |
---|
1236 | info_block['max'] = maxval |
---|
1237 | maxloc = Q.get_maximum_location() |
---|
1238 | info_block['max_location'] = maxloc |
---|
1239 | info_block['max_time'] = self.time |
---|
1240 | |
---|
1241 | # Update minimum |
---|
1242 | minval = Q.get_minimum_value(self.monitor_indices) |
---|
1243 | if info_block['min'] is None or \ |
---|
1244 | minval < info_block['min'] - epsilon: |
---|
1245 | info_block['min'] = minval |
---|
1246 | minloc = Q.get_minimum_location() |
---|
1247 | info_block['min_location'] = minloc |
---|
1248 | info_block['min_time'] = self.time |
---|
1249 | |
---|
1250 | ## |
---|
1251 | # @brief Return string with statistics about quantities |
---|
1252 | # @param precision A format string to use for float values. |
---|
1253 | # @return The stats string. |
---|
1254 | def quantity_statistics(self, precision='%.4f'): |
---|
1255 | """Return string with statistics about quantities for |
---|
1256 | printing or logging |
---|
1257 | |
---|
1258 | Quantities reported are specified through method |
---|
1259 | |
---|
1260 | set_quantities_to_be_monitored |
---|
1261 | """ |
---|
1262 | |
---|
1263 | maxlen = 128 # Max length of polygon string representation |
---|
1264 | |
---|
1265 | # Output statistics |
---|
1266 | msg = 'Monitored quantities at time %.4f:\n' % self.get_time() |
---|
1267 | if self.monitor_polygon is not None: |
---|
1268 | p_str = str(self.monitor_polygon) |
---|
1269 | msg += '- Restricted by polygon: %s' % p_str[:maxlen] |
---|
1270 | if len(p_str) >= maxlen: |
---|
1271 | msg += '...\n' |
---|
1272 | else: |
---|
1273 | msg += '\n' |
---|
1274 | |
---|
1275 | if self.monitor_time_interval is not None: |
---|
1276 | msg += '- Restricted by time interval: %s\n' \ |
---|
1277 | % str(self.monitor_time_interval) |
---|
1278 | time_interval_start = self.monitor_time_interval[0] |
---|
1279 | else: |
---|
1280 | time_interval_start = 0.0 |
---|
1281 | |
---|
1282 | for quantity_name, info in self.quantities_to_be_monitored.items(): |
---|
1283 | msg += ' %s:\n' % quantity_name |
---|
1284 | |
---|
1285 | msg += ' values since time = %.2f in [%s, %s]\n' \ |
---|
1286 | % (time_interval_start, |
---|
1287 | get_textual_float(info['min'], precision), |
---|
1288 | get_textual_float(info['max'], precision)) |
---|
1289 | |
---|
1290 | msg += ' minimum attained at time = %s, location = %s\n' \ |
---|
1291 | % (get_textual_float(info['min_time'], precision), |
---|
1292 | get_textual_float(info['min_location'], precision)) |
---|
1293 | |
---|
1294 | msg += ' maximum attained at time = %s, location = %s\n' \ |
---|
1295 | % (get_textual_float(info['max_time'], precision), |
---|
1296 | get_textual_float(info['max_location'], precision)) |
---|
1297 | |
---|
1298 | return msg |
---|
1299 | |
---|
1300 | ## |
---|
1301 | # @brief Get the timestep method. |
---|
1302 | # @return The timestep method. One of 'euler', 'rk2' or 'rk3' or 1, 2, 3. |
---|
1303 | def get_timestepping_method(self): |
---|
1304 | return self.timestepping_method |
---|
1305 | |
---|
1306 | ## |
---|
1307 | # @brief Set the tmestep method to be used. |
---|
1308 | # @param timestepping_method One of 'euler', 'rk2' or 'rk3'. |
---|
1309 | # @note Raises exception of method not known. |
---|
1310 | def set_timestepping_method(self, timestepping_method): |
---|
1311 | methods = ['euler', 'rk2', 'rk3'] |
---|
1312 | if timestepping_method in methods: |
---|
1313 | self.timestepping_method = timestepping_method |
---|
1314 | return |
---|
1315 | if timestepping_method in [1,2,3]: |
---|
1316 | self.timetepping_method = methods[timestepping_method-1] |
---|
1317 | return |
---|
1318 | |
---|
1319 | msg = '%s is an incorrect timestepping type' % timestepping_method |
---|
1320 | raise Exception, msg |
---|
1321 | |
---|
1322 | ## |
---|
1323 | # @brief Get the Domain simulation name. |
---|
1324 | # @return The simulation name string. |
---|
1325 | def get_name(self): |
---|
1326 | return self.simulation_name |
---|
1327 | |
---|
1328 | ## |
---|
1329 | # @brief Set the simulation name. |
---|
1330 | # @param name The name of the simulation. |
---|
1331 | # @note The simulation name is also used for the output .sww file. |
---|
1332 | def set_name(self, name): |
---|
1333 | """Assign a name to this simulation. |
---|
1334 | This will be used to identify the output sww file. |
---|
1335 | """ |
---|
1336 | |
---|
1337 | # remove any '.sww' end |
---|
1338 | if name.endswith('.sww'): |
---|
1339 | name = name[:-4] |
---|
1340 | |
---|
1341 | self.simulation_name = name |
---|
1342 | |
---|
1343 | ## |
---|
1344 | # @brief Get data directory path. |
---|
1345 | # @return The data directory path string. |
---|
1346 | def get_datadir(self): |
---|
1347 | return self.datadir |
---|
1348 | |
---|
1349 | ## |
---|
1350 | # @brief Set data directory path. |
---|
1351 | # @param name The data directory path string. |
---|
1352 | def set_datadir(self, name): |
---|
1353 | self.datadir = name |
---|
1354 | |
---|
1355 | ## |
---|
1356 | # @brief Get the start time value. |
---|
1357 | # @return The start time value (float). |
---|
1358 | def get_starttime(self): |
---|
1359 | return self.starttime |
---|
1360 | |
---|
1361 | ## |
---|
1362 | # @brief Set the start time value. |
---|
1363 | # @param time The start time value. |
---|
1364 | def set_starttime(self, time): |
---|
1365 | self.starttime = float(time) |
---|
1366 | |
---|
1367 | ################################################################################ |
---|
1368 | # Main components of evolve |
---|
1369 | ################################################################################ |
---|
1370 | |
---|
1371 | ## |
---|
1372 | # @brief Evolve the model through time. |
---|
1373 | # @param yieldstep Interval between yields where results are stored, etc. |
---|
1374 | # @param finaltime Time where simulation should end. |
---|
1375 | # @param duration Duration of simulation. |
---|
1376 | # @param skip_initial_step If True, skip the first yield step. |
---|
1377 | def evolve(self, yieldstep=None, |
---|
1378 | finaltime=None, |
---|
1379 | duration=None, |
---|
1380 | skip_initial_step=False): |
---|
1381 | """Evolve model through time starting from self.starttime. |
---|
1382 | |
---|
1383 | yieldstep: Interval between yields where results are stored, |
---|
1384 | statistics written and domain inspected or |
---|
1385 | possibly modified. If omitted the internal predefined |
---|
1386 | max timestep is used. |
---|
1387 | Internally, smaller timesteps may be taken. |
---|
1388 | |
---|
1389 | duration: Duration of simulation |
---|
1390 | |
---|
1391 | finaltime: Time where simulation should end. This is currently |
---|
1392 | relative time. So it's the same as duration. |
---|
1393 | |
---|
1394 | If both duration and finaltime are given an exception is thrown. |
---|
1395 | |
---|
1396 | skip_initial_step: Boolean flag that decides whether the first |
---|
1397 | yield step is skipped or not. This is useful for example to avoid |
---|
1398 | duplicate steps when multiple evolve processes are dove tailed. |
---|
1399 | |
---|
1400 | Evolve is implemented as a generator and is to be called as such, e.g. |
---|
1401 | |
---|
1402 | for t in domain.evolve(yieldstep, finaltime): |
---|
1403 | <Do something with domain and t> |
---|
1404 | |
---|
1405 | All times are given in seconds |
---|
1406 | """ |
---|
1407 | |
---|
1408 | from anuga.config import epsilon |
---|
1409 | |
---|
1410 | # FIXME: Maybe lump into a larger check prior to evolving |
---|
1411 | msg = ('Boundary tags must be bound to boundary objects before ' |
---|
1412 | 'evolving system, ' |
---|
1413 | 'e.g. using the method set_boundary.\n' |
---|
1414 | 'This system has the boundary tags %s ' |
---|
1415 | % self.get_boundary_tags()) |
---|
1416 | assert hasattr(self, 'boundary_objects'), msg |
---|
1417 | |
---|
1418 | if yieldstep is None: |
---|
1419 | yieldstep = self.evolve_max_timestep |
---|
1420 | else: |
---|
1421 | yieldstep = float(yieldstep) |
---|
1422 | |
---|
1423 | self._order_ = self.default_order |
---|
1424 | |
---|
1425 | if finaltime is not None and duration is not None: |
---|
1426 | msg = 'Only one of finaltime and duration may be specified' |
---|
1427 | raise Exception, msg |
---|
1428 | else: |
---|
1429 | if finaltime is not None: |
---|
1430 | self.finaltime = float(finaltime) |
---|
1431 | if duration is not None: |
---|
1432 | self.finaltime = self.starttime + float(duration) |
---|
1433 | |
---|
1434 | N = len(self) # Number of triangles |
---|
1435 | self.yieldtime = self.time + yieldstep # set next yield time |
---|
1436 | |
---|
1437 | # Initialise interval of timestep sizes (for reporting only) |
---|
1438 | # Note that we set recorded_min_timestep to be large so that it comes |
---|
1439 | # down through the evolution, similarly recorded_max_timestep |
---|
1440 | self.recorded_min_timestep = self.evolve_max_timestep |
---|
1441 | self.recorded_max_timestep = self.evolve_min_timestep |
---|
1442 | self.number_of_steps = 0 |
---|
1443 | self.number_of_first_order_steps = 0 |
---|
1444 | |
---|
1445 | # Update ghosts |
---|
1446 | self.update_ghosts() |
---|
1447 | |
---|
1448 | # Initial update of vertex and edge values |
---|
1449 | self.distribute_to_vertices_and_edges() |
---|
1450 | |
---|
1451 | # Initial update boundary values |
---|
1452 | self.update_boundary() |
---|
1453 | |
---|
1454 | # Update extrema if necessary (for reporting) |
---|
1455 | self.update_extrema() |
---|
1456 | |
---|
1457 | |
---|
1458 | |
---|
1459 | # Or maybe restore from latest checkpoint |
---|
1460 | if self.checkpoint is True: |
---|
1461 | self.goto_latest_checkpoint() |
---|
1462 | |
---|
1463 | if skip_initial_step is False: |
---|
1464 | yield(self.time) # Yield initial values |
---|
1465 | |
---|
1466 | while True: |
---|
1467 | |
---|
1468 | initial_time = self.time |
---|
1469 | |
---|
1470 | #========================================== |
---|
1471 | # Apply fluid flow fractional step |
---|
1472 | #========================================== |
---|
1473 | if self.get_timestepping_method() == 'euler': |
---|
1474 | self.evolve_one_euler_step(yieldstep, finaltime) |
---|
1475 | |
---|
1476 | elif self.get_timestepping_method() == 'rk2': |
---|
1477 | self.evolve_one_rk2_step(yieldstep, finaltime) |
---|
1478 | |
---|
1479 | elif self.get_timestepping_method() == 'rk3': |
---|
1480 | self.evolve_one_rk3_step(yieldstep, finaltime) |
---|
1481 | |
---|
1482 | |
---|
1483 | #========================================== |
---|
1484 | # Apply other fractional steps |
---|
1485 | #========================================== |
---|
1486 | #self.apply_fractional_steps() |
---|
1487 | |
---|
1488 | #========================================== |
---|
1489 | # Centroid Values of variables should be ok, |
---|
1490 | # so now setup quantites etc for output |
---|
1491 | #========================================== |
---|
1492 | |
---|
1493 | # Update time |
---|
1494 | self.time = initial_time + self.timestep |
---|
1495 | |
---|
1496 | # Update vertex and edge values |
---|
1497 | self.distribute_to_vertices_and_edges() |
---|
1498 | |
---|
1499 | # Update boundary values |
---|
1500 | self.update_boundary() |
---|
1501 | |
---|
1502 | # Update extrema if necessary (for reporting) |
---|
1503 | self.update_extrema() |
---|
1504 | |
---|
1505 | self.number_of_steps += 1 |
---|
1506 | if self._order_ == 1: |
---|
1507 | self.number_of_first_order_steps += 1 |
---|
1508 | |
---|
1509 | # Yield results |
---|
1510 | if finaltime is not None and self.time >= finaltime-epsilon: |
---|
1511 | if self.time > finaltime: |
---|
1512 | # FIXME (Ole, 30 April 2006): Do we need this check? |
---|
1513 | # Probably not (Ole, 18 September 2008). |
---|
1514 | # Now changed to Exception. |
---|
1515 | msg = ('WARNING (domain.py): time overshot finaltime. ' |
---|
1516 | 'Contact Ole.Nielsen@ga.gov.au') |
---|
1517 | raise Exception, msg |
---|
1518 | |
---|
1519 | # Yield final time and stop |
---|
1520 | self.time = finaltime |
---|
1521 | yield(self.time) |
---|
1522 | break |
---|
1523 | |
---|
1524 | # if we are at the next yield point |
---|
1525 | if self.time >= self.yieldtime: |
---|
1526 | # Yield (intermediate) time and allow inspection of domain |
---|
1527 | if self.checkpoint is True: |
---|
1528 | self.store_checkpoint() |
---|
1529 | self.delete_old_checkpoints() |
---|
1530 | |
---|
1531 | # Pass control on to outer loop for more specific actions |
---|
1532 | yield(self.time) |
---|
1533 | |
---|
1534 | # Reinitialise |
---|
1535 | self.yieldtime += yieldstep # move to next yield |
---|
1536 | self.recorded_min_timestep = self.evolve_max_timestep |
---|
1537 | self.recorded_max_timestep = self.evolve_min_timestep |
---|
1538 | self.number_of_steps = 0 |
---|
1539 | self.number_of_first_order_steps = 0 |
---|
1540 | self.max_speed = num.zeros(N, num.float) |
---|
1541 | |
---|
1542 | ## |
---|
1543 | # @brief 'Euler' time step method. |
---|
1544 | # @param yieldstep The reporting time step. |
---|
1545 | # @param finaltime The simulation final time. |
---|
1546 | def evolve_one_euler_step(self, yieldstep, finaltime): |
---|
1547 | """One Euler Time Step |
---|
1548 | Q^{n+1} = E(h) Q^n |
---|
1549 | |
---|
1550 | Assumes that centroid values have been extrapolated to vertices and edges |
---|
1551 | """ |
---|
1552 | |
---|
1553 | # Compute fluxes across each element edge |
---|
1554 | self.compute_fluxes() |
---|
1555 | |
---|
1556 | # Compute forcing terms |
---|
1557 | self.compute_forcing_terms() |
---|
1558 | |
---|
1559 | # Update timestep to fit yieldstep and finaltime |
---|
1560 | self.update_timestep(yieldstep, finaltime) |
---|
1561 | |
---|
1562 | # Update conserved quantities |
---|
1563 | self.update_conserved_quantities() |
---|
1564 | |
---|
1565 | # Update ghosts |
---|
1566 | self.update_ghosts() |
---|
1567 | |
---|
1568 | |
---|
1569 | |
---|
1570 | |
---|
1571 | ## |
---|
1572 | # @brief 'rk2' time step method. |
---|
1573 | # @param yieldstep The reporting time step. |
---|
1574 | # @param finaltime The simulation final time. |
---|
1575 | def evolve_one_rk2_step(self, yieldstep, finaltime): |
---|
1576 | """One 2nd order RK timestep |
---|
1577 | Q^{n+1} = 0.5 Q^n + 0.5 E(h)^2 Q^n |
---|
1578 | """ |
---|
1579 | |
---|
1580 | # Save initial initial conserved quantities values |
---|
1581 | self.backup_conserved_quantities() |
---|
1582 | |
---|
1583 | ###### |
---|
1584 | # First euler step |
---|
1585 | ###### |
---|
1586 | |
---|
1587 | # Compute fluxes across each element edge |
---|
1588 | self.compute_fluxes() |
---|
1589 | |
---|
1590 | # Compute forcing terms |
---|
1591 | self.compute_forcing_terms() |
---|
1592 | |
---|
1593 | # Update timestep to fit yieldstep and finaltime |
---|
1594 | self.update_timestep(yieldstep, finaltime) |
---|
1595 | |
---|
1596 | # Update conserved quantities |
---|
1597 | self.update_conserved_quantities() |
---|
1598 | |
---|
1599 | # Update ghosts |
---|
1600 | self.update_ghosts() |
---|
1601 | |
---|
1602 | # Update time |
---|
1603 | self.time += self.timestep |
---|
1604 | |
---|
1605 | # Update vertex and edge values |
---|
1606 | self.distribute_to_vertices_and_edges() |
---|
1607 | |
---|
1608 | # Update boundary values |
---|
1609 | self.update_boundary() |
---|
1610 | |
---|
1611 | ###### |
---|
1612 | # Second Euler step using the same timestep |
---|
1613 | # calculated in the first step. Might lead to |
---|
1614 | # stability problems but we have not seen any |
---|
1615 | # example. |
---|
1616 | ###### |
---|
1617 | |
---|
1618 | # Compute fluxes across each element edge |
---|
1619 | self.compute_fluxes() |
---|
1620 | |
---|
1621 | # Compute forcing terms |
---|
1622 | self.compute_forcing_terms() |
---|
1623 | |
---|
1624 | # Update conserved quantities |
---|
1625 | self.update_conserved_quantities() |
---|
1626 | |
---|
1627 | ###### |
---|
1628 | # Combine initial and final values |
---|
1629 | # of conserved quantities and cleanup |
---|
1630 | ###### |
---|
1631 | |
---|
1632 | # Combine steps |
---|
1633 | self.saxpy_conserved_quantities(0.5, 0.5) |
---|
1634 | |
---|
1635 | # Update ghosts |
---|
1636 | self.update_ghosts() |
---|
1637 | |
---|
1638 | |
---|
1639 | ## |
---|
1640 | # @brief 'rk3' time step method. |
---|
1641 | # @param yieldstep The reporting time step. |
---|
1642 | # @param finaltime The simulation final time. |
---|
1643 | def evolve_one_rk3_step(self, yieldstep, finaltime): |
---|
1644 | """One 3rd order RK timestep |
---|
1645 | Q^(1) = 3/4 Q^n + 1/4 E(h)^2 Q^n (at time t^n + h/2) |
---|
1646 | Q^{n+1} = 1/3 Q^n + 2/3 E(h) Q^(1) (at time t^{n+1}) |
---|
1647 | """ |
---|
1648 | |
---|
1649 | # Save initial initial conserved quantities values |
---|
1650 | self.backup_conserved_quantities() |
---|
1651 | |
---|
1652 | initial_time = self.time |
---|
1653 | |
---|
1654 | ###### |
---|
1655 | # First euler step |
---|
1656 | ###### |
---|
1657 | |
---|
1658 | # Compute fluxes across each element edge |
---|
1659 | self.compute_fluxes() |
---|
1660 | |
---|
1661 | # Compute forcing terms |
---|
1662 | self.compute_forcing_terms() |
---|
1663 | |
---|
1664 | # Update timestep to fit yieldstep and finaltime |
---|
1665 | self.update_timestep(yieldstep, finaltime) |
---|
1666 | |
---|
1667 | # Update conserved quantities |
---|
1668 | self.update_conserved_quantities() |
---|
1669 | |
---|
1670 | # Update ghosts |
---|
1671 | self.update_ghosts() |
---|
1672 | |
---|
1673 | # Update time |
---|
1674 | self.time += self.timestep |
---|
1675 | |
---|
1676 | # Update vertex and edge values |
---|
1677 | self.distribute_to_vertices_and_edges() |
---|
1678 | |
---|
1679 | # Update boundary values |
---|
1680 | self.update_boundary() |
---|
1681 | |
---|
1682 | ###### |
---|
1683 | # Second Euler step using the same timestep |
---|
1684 | # calculated in the first step. Might lead to |
---|
1685 | # stability problems but we have not seen any |
---|
1686 | # example. |
---|
1687 | ###### |
---|
1688 | |
---|
1689 | # Compute fluxes across each element edge |
---|
1690 | self.compute_fluxes() |
---|
1691 | |
---|
1692 | # Compute forcing terms |
---|
1693 | self.compute_forcing_terms() |
---|
1694 | |
---|
1695 | # Update conserved quantities |
---|
1696 | self.update_conserved_quantities() |
---|
1697 | |
---|
1698 | ###### |
---|
1699 | # Combine steps to obtain intermediate |
---|
1700 | # solution at time t^n + 0.5 h |
---|
1701 | ###### |
---|
1702 | |
---|
1703 | # Combine steps |
---|
1704 | self.saxpy_conserved_quantities(0.25, 0.75) |
---|
1705 | |
---|
1706 | # Update ghosts |
---|
1707 | self.update_ghosts() |
---|
1708 | |
---|
1709 | # Set substep time |
---|
1710 | self.time = initial_time + self.timestep*0.5 |
---|
1711 | |
---|
1712 | # Update vertex and edge values |
---|
1713 | self.distribute_to_vertices_and_edges() |
---|
1714 | |
---|
1715 | # Update boundary values |
---|
1716 | self.update_boundary() |
---|
1717 | |
---|
1718 | ###### |
---|
1719 | # Third Euler step |
---|
1720 | ###### |
---|
1721 | |
---|
1722 | # Compute fluxes across each element edge |
---|
1723 | self.compute_fluxes() |
---|
1724 | |
---|
1725 | # Compute forcing terms |
---|
1726 | self.compute_forcing_terms() |
---|
1727 | |
---|
1728 | # Update conserved quantities |
---|
1729 | self.update_conserved_quantities() |
---|
1730 | |
---|
1731 | ###### |
---|
1732 | # Combine final and initial values |
---|
1733 | # and cleanup |
---|
1734 | ###### |
---|
1735 | |
---|
1736 | # Combine steps |
---|
1737 | self.saxpy_conserved_quantities(2.0/3.0, 1.0/3.0) |
---|
1738 | |
---|
1739 | # Update ghosts |
---|
1740 | self.update_ghosts() |
---|
1741 | |
---|
1742 | # Set new time |
---|
1743 | self.time = initial_time + self.timestep |
---|
1744 | |
---|
1745 | |
---|
1746 | ## |
---|
1747 | # @brief Evolve simulation to a final time. |
---|
1748 | # @param finaltime Sinulation final time. |
---|
1749 | def evolve_to_end(self, finaltime=1.0): |
---|
1750 | """Iterate evolve all the way to the end.""" |
---|
1751 | |
---|
1752 | for _ in self.evolve(yieldstep=None, finaltime=finaltime): |
---|
1753 | pass |
---|
1754 | |
---|
1755 | ## |
---|
1756 | # @brief Backup conserved quantities |
---|
1757 | def backup_conserved_quantities(self): |
---|
1758 | |
---|
1759 | # Backup conserved_quantities centroid values |
---|
1760 | for name in self.conserved_quantities: |
---|
1761 | Q = self.quantities[name] |
---|
1762 | Q.backup_centroid_values() |
---|
1763 | |
---|
1764 | ## |
---|
1765 | # @brief Combines current C and saved centroid values S as C = aC + bS |
---|
1766 | # @param a factor in combination |
---|
1767 | # @param b factor in combination |
---|
1768 | def saxpy_conserved_quantities(self, a, b): |
---|
1769 | |
---|
1770 | # Backup conserved_quantities centroid values |
---|
1771 | for name in self.conserved_quantities: |
---|
1772 | Q = self.quantities[name] |
---|
1773 | Q.saxpy_centroid_values(a, b) |
---|
1774 | |
---|
1775 | |
---|
1776 | |
---|
1777 | |
---|
1778 | ## |
---|
1779 | # @brief Mapping between conserved quantites and evolved quantities |
---|
1780 | # @param Input: q_cons array of conserved quantity values |
---|
1781 | # @param Input: q_evol array of current evolved quantity values |
---|
1782 | # @note Output: Updated q_evol array |
---|
1783 | def conserved_values_to_evolved_values(self, q_cons, q_evol): |
---|
1784 | """Needs to be overridden by Domain subclass |
---|
1785 | """ |
---|
1786 | |
---|
1787 | if len(q_cons) == len(q_evol): |
---|
1788 | q_evol[:] = q_cons |
---|
1789 | else: |
---|
1790 | msg = 'Method conserved_values_to_evolved_values must be overridden' |
---|
1791 | msg += ' by Domain subclass' |
---|
1792 | raise Exception, msg |
---|
1793 | |
---|
1794 | return q_evol |
---|
1795 | |
---|
1796 | ## |
---|
1797 | # @brief Update boundary values for all conserved quantities. |
---|
1798 | def update_boundary(self): |
---|
1799 | """Go through list of boundary objects and update boundary values |
---|
1800 | for all conserved quantities on boundary. |
---|
1801 | It is assumed that the ordering of conserved quantities is |
---|
1802 | consistent between the domain and the boundary object, i.e. |
---|
1803 | the jth element of vector q must correspond to the jth conserved |
---|
1804 | quantity in domain. |
---|
1805 | """ |
---|
1806 | |
---|
1807 | # FIXME: Update only those that change (if that can be worked out) |
---|
1808 | # FIXME: Boundary objects should not include ghost nodes. |
---|
1809 | for i, ((vol_id, edge_id), B) in enumerate(self.boundary_objects): |
---|
1810 | if B is None: |
---|
1811 | log.critical('WARNING: Ignored boundary segment (None)') |
---|
1812 | else: |
---|
1813 | q_bdry = B.evaluate(vol_id, edge_id) |
---|
1814 | |
---|
1815 | if len(q_bdry) == len(self.evolved_quantities): |
---|
1816 | # conserved and evolved quantities are the same |
---|
1817 | q_evol = q_bdry |
---|
1818 | elif len(q_bdry) == len(self.conserved_quantities): |
---|
1819 | # boundary just returns conserved quantities |
---|
1820 | # Need to calculate all the evolved quantities |
---|
1821 | # Use default conversion |
---|
1822 | |
---|
1823 | q_evol = self.get_evolved_quantities(vol_id, edge = edge_id) |
---|
1824 | |
---|
1825 | q_evol = self.conserved_values_to_evolved_values \ |
---|
1826 | (q_bdry, q_evol) |
---|
1827 | else: |
---|
1828 | msg = 'Boundary must return array of either conserved' |
---|
1829 | msg += ' or evolved quantities' |
---|
1830 | raise Exception, msg |
---|
1831 | |
---|
1832 | for j, name in enumerate(self.evolved_quantities): |
---|
1833 | Q = self.quantities[name] |
---|
1834 | Q.boundary_values[i] = q_evol[j] |
---|
1835 | |
---|
1836 | ## |
---|
1837 | # @brief Compute fluxes. |
---|
1838 | # @note MUST BE OVERRIDEN IN SUBCLASS! |
---|
1839 | def compute_fluxes(self): |
---|
1840 | msg = 'Method compute_fluxes must be overridden by Domain subclass' |
---|
1841 | raise Exception, msg |
---|
1842 | |
---|
1843 | ## |
---|
1844 | # @brief |
---|
1845 | # @param yieldstep |
---|
1846 | # @param finaltime |
---|
1847 | def update_timestep(self, yieldstep, finaltime): |
---|
1848 | |
---|
1849 | # Protect against degenerate timesteps arising from isolated |
---|
1850 | # triangles |
---|
1851 | self.apply_protection_against_isolated_degenerate_timesteps() |
---|
1852 | |
---|
1853 | # self.timestep is calculated from speed of characteristics |
---|
1854 | # Apply CFL condition here |
---|
1855 | timestep = min(self.CFL*self.flux_timestep, self.evolve_max_timestep) |
---|
1856 | |
---|
1857 | # Record maximal and minimal values of timestep for reporting |
---|
1858 | self.recorded_max_timestep = max(timestep, self.recorded_max_timestep) |
---|
1859 | self.recorded_min_timestep = min(timestep, self.recorded_min_timestep) |
---|
1860 | |
---|
1861 | # Protect against degenerate time steps |
---|
1862 | if timestep < self.evolve_min_timestep: |
---|
1863 | # Number of consecutive small steps taken b4 taking action |
---|
1864 | self.smallsteps += 1 |
---|
1865 | |
---|
1866 | if self.smallsteps > self.max_smallsteps: |
---|
1867 | self.smallsteps = 0 # Reset |
---|
1868 | |
---|
1869 | if self._order_ == 1: |
---|
1870 | msg = 'WARNING: Too small timestep %.16f reached ' \ |
---|
1871 | % timestep |
---|
1872 | msg += 'even after %d steps of 1 order scheme' \ |
---|
1873 | % self.max_smallsteps |
---|
1874 | log.critical(msg) |
---|
1875 | timestep = self.evolve_min_timestep # Try enforce min_step |
---|
1876 | |
---|
1877 | stats = self.timestepping_statistics(track_speeds=True) |
---|
1878 | log.critical(stats) |
---|
1879 | |
---|
1880 | raise Exception, msg |
---|
1881 | else: |
---|
1882 | # Try to overcome situation by switching to 1 order |
---|
1883 | self._order_ = 1 |
---|
1884 | else: |
---|
1885 | self.smallsteps = 0 |
---|
1886 | if self._order_ == 1 and self.default_order == 2: |
---|
1887 | self._order_ = 2 |
---|
1888 | |
---|
1889 | # Ensure that final time is not exceeded |
---|
1890 | if finaltime is not None and self.time + timestep > finaltime : |
---|
1891 | timestep = finaltime-self.time |
---|
1892 | |
---|
1893 | # Ensure that model time is aligned with yieldsteps |
---|
1894 | if self.time + timestep > self.yieldtime: |
---|
1895 | timestep = self.yieldtime - self.time |
---|
1896 | |
---|
1897 | self.timestep = timestep |
---|
1898 | |
---|
1899 | ## |
---|
1900 | # @brief Compute forcing terms, if any. |
---|
1901 | def compute_forcing_terms(self): |
---|
1902 | """If there are any forcing functions driving the system |
---|
1903 | they should be defined in Domain subclass and appended to |
---|
1904 | the list self.forcing_terms |
---|
1905 | """ |
---|
1906 | |
---|
1907 | # The parameter self.flux_timestep should be updated |
---|
1908 | # by the forcing_terms to ensure stability |
---|
1909 | |
---|
1910 | for f in self.forcing_terms: |
---|
1911 | f(self) |
---|
1912 | |
---|
1913 | |
---|
1914 | ## |
---|
1915 | # @brief Update vectors of conserved quantities. |
---|
1916 | def update_conserved_quantities(self): |
---|
1917 | """Update vectors of conserved quantities using previously |
---|
1918 | computed fluxes and specified forcing functions. |
---|
1919 | """ |
---|
1920 | |
---|
1921 | N = len(self) # Number_of_triangles |
---|
1922 | d = len(self.conserved_quantities) |
---|
1923 | |
---|
1924 | timestep = self.timestep |
---|
1925 | |
---|
1926 | |
---|
1927 | # Update conserved_quantities |
---|
1928 | for name in self.conserved_quantities: |
---|
1929 | Q = self.quantities[name] |
---|
1930 | Q.update(timestep) |
---|
1931 | |
---|
1932 | # Note that Q.explicit_update is reset by compute_fluxes |
---|
1933 | # Where is Q.semi_implicit_update reset? |
---|
1934 | # It is reset in quantity_ext.c |
---|
1935 | |
---|
1936 | ## |
---|
1937 | # @brief Sequential update of ghost cells |
---|
1938 | def update_ghosts(self): |
---|
1939 | # We must send the information from the full cells and |
---|
1940 | # receive the information for the ghost cells |
---|
1941 | # We have a list with ghosts expecting updates |
---|
1942 | |
---|
1943 | #Update of ghost cells |
---|
1944 | iproc = self.processor |
---|
1945 | if self.full_send_dict.has_key(iproc): |
---|
1946 | |
---|
1947 | # now store full as local id, global id, value |
---|
1948 | Idf = self.full_send_dict[iproc][0] |
---|
1949 | |
---|
1950 | # now store ghost as local id, global id, value |
---|
1951 | Idg = self.ghost_recv_dict[iproc][0] |
---|
1952 | |
---|
1953 | for i, q in enumerate(self.conserved_quantities): |
---|
1954 | Q_cv = self.quantities[q].centroid_values |
---|
1955 | num.put(Q_cv, Idg, num.take(Q_cv, Idf, axis=0)) |
---|
1956 | |
---|
1957 | |
---|
1958 | ## |
---|
1959 | # @brief Extrapolate conserved quantities from centroid to vertices |
---|
1960 | # and edge-midpoints for each volume. |
---|
1961 | def distribute_to_vertices_and_edges(self): |
---|
1962 | """Extrapolate conserved quantities from centroid to |
---|
1963 | vertices and edge-midpoints for each volume |
---|
1964 | |
---|
1965 | Default implementation is straight first order, |
---|
1966 | i.e. constant values throughout each element and |
---|
1967 | no reference to non-conserved quantities. |
---|
1968 | """ |
---|
1969 | |
---|
1970 | for name in self.conserved_quantities: |
---|
1971 | Q = self.quantities[name] |
---|
1972 | if self._order_ == 1: |
---|
1973 | Q.extrapolate_first_order() |
---|
1974 | elif self._order_ == 2: |
---|
1975 | Q.extrapolate_second_order() |
---|
1976 | else: |
---|
1977 | raise Exception, 'Unknown order' |
---|
1978 | |
---|
1979 | ## |
---|
1980 | # @brief Calculate the norm of the centroid values of a specific quantity, |
---|
1981 | # using normfunc. |
---|
1982 | # @param quantity |
---|
1983 | # @param normfunc |
---|
1984 | def centroid_norm(self, quantity, normfunc): |
---|
1985 | """Calculate the norm of the centroid values of a specific quantity, |
---|
1986 | using normfunc. |
---|
1987 | |
---|
1988 | normfunc should take a list to a float. |
---|
1989 | |
---|
1990 | common normfuncs are provided in the module utilities.norms |
---|
1991 | """ |
---|
1992 | |
---|
1993 | return normfunc(self.quantities[quantity].centroid_values) |
---|
1994 | |
---|
1995 | |
---|
1996 | |
---|
1997 | def apply_protection_against_isolated_degenerate_timesteps(self): |
---|
1998 | |
---|
1999 | # FIXME (Steve): This should be in shallow_water as it assumes x and y |
---|
2000 | # momentum |
---|
2001 | if self.protect_against_isolated_degenerate_timesteps is False: |
---|
2002 | return |
---|
2003 | |
---|
2004 | # FIXME (Ole): Make this configurable |
---|
2005 | if num.max(self.max_speed) < 10.0: |
---|
2006 | return |
---|
2007 | |
---|
2008 | # Setup 10 bins for speed histogram |
---|
2009 | from anuga.utilities.numerical_tools import histogram, create_bins |
---|
2010 | |
---|
2011 | bins = create_bins(self.max_speed, 10) |
---|
2012 | hist = histogram(self.max_speed, bins) |
---|
2013 | |
---|
2014 | # Look for characteristic signature |
---|
2015 | if len(hist) > 1 and hist[-1] > 0 and \ |
---|
2016 | hist[4] == hist[5] == hist[6] == hist[7] == hist[8] == 0: |
---|
2017 | # Danger of isolated degenerate triangles |
---|
2018 | |
---|
2019 | # Find triangles in last bin |
---|
2020 | # FIXME - speed up using numeric package |
---|
2021 | d = 0 |
---|
2022 | for i in range(self.number_of_full_triangles): |
---|
2023 | if self.max_speed[i] > bins[-1]: |
---|
2024 | msg = 'Time=%f: Ignoring isolated high ' % self.time |
---|
2025 | msg += 'speed triangle ' |
---|
2026 | msg += '#%d of %d with max speed=%f' \ |
---|
2027 | % (i, self.number_of_full_triangles, self.max_speed[i]) |
---|
2028 | |
---|
2029 | self.get_quantity('xmomentum').\ |
---|
2030 | set_values(0.0, indices=[i]) |
---|
2031 | self.get_quantity('ymomentum').\ |
---|
2032 | set_values(0.0, indices=[i]) |
---|
2033 | self.max_speed[i]=0.0 |
---|
2034 | d += 1 |
---|
2035 | |
---|
2036 | |
---|
2037 | ###### |
---|
2038 | # Initialise module |
---|
2039 | ###### |
---|
2040 | |
---|
2041 | # Optimisation with psyco |
---|
2042 | #from anuga.config import use_psyco |
---|
2043 | |
---|
2044 | #if use_psyco: |
---|
2045 | #try: |
---|
2046 | #import psyco |
---|
2047 | #except: |
---|
2048 | #import os |
---|
2049 | #if os.name == 'posix' and os.uname()[4] in ['x86_64', 'ia64']: |
---|
2050 | #pass |
---|
2051 | ## Psyco isn't supported on 64 bit systems, but it doesn't matter |
---|
2052 | #else: |
---|
2053 | #log.critical('WARNING: psyco (speedup) could not be imported, ' |
---|
2054 | #'you may want to consider installing it') |
---|
2055 | #else: |
---|
2056 | #psyco.bind(Generic_Domain.update_boundary) |
---|
2057 | ##psyco.bind(Domain.update_timestep) # Not worth it |
---|
2058 | #psyco.bind(Generic_Domain.update_conserved_quantities) |
---|
2059 | #psyco.bind(Generic_Domain.distribute_to_vertices_and_edges) |
---|
2060 | |
---|
2061 | |
---|
2062 | if __name__ == "__main__": |
---|
2063 | pass |
---|