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