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