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 | from Numeric import allclose, argmax |
---|
11 | from anuga.config import epsilon |
---|
12 | |
---|
13 | from anuga.abstract_2d_finite_volumes.neighbour_mesh import Mesh |
---|
14 | from anuga.abstract_2d_finite_volumes.generic_boundary_conditions\ |
---|
15 | import Boundary |
---|
16 | from anuga.abstract_2d_finite_volumes.generic_boundary_conditions\ |
---|
17 | import File_boundary |
---|
18 | from anuga.abstract_2d_finite_volumes.generic_boundary_conditions\ |
---|
19 | import Dirichlet_boundary |
---|
20 | from anuga.abstract_2d_finite_volumes.generic_boundary_conditions\ |
---|
21 | import Time_boundary |
---|
22 | from anuga.abstract_2d_finite_volumes.generic_boundary_conditions\ |
---|
23 | import Transmissive_boundary |
---|
24 | |
---|
25 | from anuga.abstract_2d_finite_volumes.pmesh2domain import pmesh_to_domain |
---|
26 | from anuga.abstract_2d_finite_volumes.region\ |
---|
27 | import Set_region as region_set_region |
---|
28 | |
---|
29 | import types |
---|
30 | |
---|
31 | class Domain(Mesh): |
---|
32 | |
---|
33 | |
---|
34 | def __init__(self, |
---|
35 | source=None, |
---|
36 | triangles=None, |
---|
37 | boundary=None, |
---|
38 | conserved_quantities=None, |
---|
39 | other_quantities=None, |
---|
40 | tagged_elements=None, |
---|
41 | geo_reference=None, |
---|
42 | use_inscribed_circle=False, |
---|
43 | mesh_filename=None, |
---|
44 | use_cache=False, |
---|
45 | verbose=False, |
---|
46 | full_send_dict=None, |
---|
47 | ghost_recv_dict=None, |
---|
48 | processor=0, |
---|
49 | numproc=1, |
---|
50 | number_of_full_nodes=None, |
---|
51 | number_of_full_triangles=None): |
---|
52 | |
---|
53 | |
---|
54 | """Instantiate generic computational Domain. |
---|
55 | |
---|
56 | Input: |
---|
57 | source: Either a mesh filename or coordinates of mesh vertices. |
---|
58 | If it is a filename values specified for triangles will |
---|
59 | be overridden. |
---|
60 | triangles: Mesh connectivity (see mesh.py for more information) |
---|
61 | boundary: See mesh.py for more information |
---|
62 | |
---|
63 | conserved_quantities: List of quantity names entering the |
---|
64 | conservation equations |
---|
65 | other_quantities: List of other quantity names |
---|
66 | |
---|
67 | tagged_elements: |
---|
68 | ... |
---|
69 | |
---|
70 | |
---|
71 | """ |
---|
72 | |
---|
73 | # Determine whether source is a mesh filename or coordinates |
---|
74 | if type(source) == types.StringType: |
---|
75 | mesh_filename = source |
---|
76 | else: |
---|
77 | coordinates = source |
---|
78 | |
---|
79 | |
---|
80 | # In case a filename has been specified, extract content |
---|
81 | if mesh_filename is not None: |
---|
82 | coordinates, triangles, boundary, vertex_quantity_dict, \ |
---|
83 | tagged_elements, geo_reference = \ |
---|
84 | pmesh_to_domain(file_name=mesh_filename, |
---|
85 | use_cache=use_cache, |
---|
86 | verbose=verbose) |
---|
87 | |
---|
88 | |
---|
89 | # Initialise underlying mesh structure |
---|
90 | Mesh.__init__(self, coordinates, triangles, |
---|
91 | boundary=boundary, |
---|
92 | tagged_elements=tagged_elements, |
---|
93 | geo_reference=geo_reference, |
---|
94 | use_inscribed_circle=use_inscribed_circle, |
---|
95 | number_of_full_nodes=number_of_full_nodes, |
---|
96 | number_of_full_triangles=number_of_full_triangles, |
---|
97 | verbose=verbose) |
---|
98 | |
---|
99 | if verbose: print 'Initialising Domain' |
---|
100 | from Numeric import zeros, Float, Int, ones |
---|
101 | from quantity import Quantity, Conserved_quantity |
---|
102 | |
---|
103 | # List of quantity names entering |
---|
104 | # the conservation equations |
---|
105 | if conserved_quantities is None: |
---|
106 | self.conserved_quantities = [] |
---|
107 | else: |
---|
108 | self.conserved_quantities = conserved_quantities |
---|
109 | |
---|
110 | # List of other quantity names |
---|
111 | if other_quantities is None: |
---|
112 | self.other_quantities = [] |
---|
113 | else: |
---|
114 | self.other_quantities = other_quantities |
---|
115 | |
---|
116 | |
---|
117 | #Build dictionary of Quantity instances keyed by quantity names |
---|
118 | self.quantities = {} |
---|
119 | |
---|
120 | #FIXME: remove later - maybe OK, though.... |
---|
121 | for name in self.conserved_quantities: |
---|
122 | self.quantities[name] = Conserved_quantity(self) |
---|
123 | for name in self.other_quantities: |
---|
124 | self.quantities[name] = Quantity(self) |
---|
125 | |
---|
126 | #Create an empty list for explicit forcing terms |
---|
127 | self.forcing_terms = [] |
---|
128 | |
---|
129 | #Setup the ghost cell communication |
---|
130 | if full_send_dict is None: |
---|
131 | self.full_send_dict = {} |
---|
132 | else: |
---|
133 | self.full_send_dict = full_send_dict |
---|
134 | |
---|
135 | # List of other quantity names |
---|
136 | if ghost_recv_dict is None: |
---|
137 | self.ghost_recv_dict = {} |
---|
138 | else: |
---|
139 | self.ghost_recv_dict = ghost_recv_dict |
---|
140 | |
---|
141 | self.processor = processor |
---|
142 | self.numproc = numproc |
---|
143 | |
---|
144 | |
---|
145 | # Setup Communication Buffers |
---|
146 | if verbose: print 'Domain: Set up communication buffers (parallel)' |
---|
147 | self.nsys = len(self.conserved_quantities) |
---|
148 | for key in self.full_send_dict: |
---|
149 | buffer_shape = self.full_send_dict[key][0].shape[0] |
---|
150 | self.full_send_dict[key].append(zeros( (buffer_shape,self.nsys) ,Float)) |
---|
151 | |
---|
152 | |
---|
153 | for key in self.ghost_recv_dict: |
---|
154 | buffer_shape = self.ghost_recv_dict[key][0].shape[0] |
---|
155 | self.ghost_recv_dict[key].append(zeros( (buffer_shape,self.nsys) ,Float)) |
---|
156 | |
---|
157 | |
---|
158 | # Setup cell full flag |
---|
159 | # =1 for full |
---|
160 | # =0 for ghost |
---|
161 | N = len(self) #number_of_elements |
---|
162 | self.tri_full_flag = ones(N, Int) |
---|
163 | for i in self.ghost_recv_dict.keys(): |
---|
164 | for id in self.ghost_recv_dict[i][0]: |
---|
165 | self.tri_full_flag[id] = 0 |
---|
166 | |
---|
167 | # Test the assumption that all full triangles are store before |
---|
168 | # the ghost triangles. |
---|
169 | assert allclose(self.tri_full_flag[:self.number_of_full_nodes],1) |
---|
170 | |
---|
171 | |
---|
172 | #Defaults |
---|
173 | from anuga.config import max_smallsteps, beta_w, beta_h, epsilon, CFL |
---|
174 | self.beta_w = beta_w |
---|
175 | self.beta_h = beta_h |
---|
176 | self.epsilon = epsilon |
---|
177 | |
---|
178 | #FIXME: Maybe have separate orders for h-limiter and w-limiter? |
---|
179 | #Or maybe get rid of order altogether and use beta_w and beta_h |
---|
180 | self.set_default_order(1) |
---|
181 | #self.default_order = 1 |
---|
182 | #self._order_ = self.default_order |
---|
183 | |
---|
184 | self.smallsteps = 0 |
---|
185 | self.max_smallsteps = max_smallsteps |
---|
186 | self.number_of_steps = 0 |
---|
187 | self.number_of_first_order_steps = 0 |
---|
188 | self.CFL = CFL |
---|
189 | |
---|
190 | self.boundary_map = None # Will be populated by set_boundary |
---|
191 | |
---|
192 | |
---|
193 | #Model time |
---|
194 | self.time = 0.0 |
---|
195 | self.finaltime = None |
---|
196 | self.min_timestep = self.max_timestep = 0.0 |
---|
197 | self.starttime = 0 #Physical starttime if any (0 is 1 Jan 1970 00:00:00) |
---|
198 | |
---|
199 | ######OBSOLETE |
---|
200 | #Origin in UTM coordinates |
---|
201 | #FIXME: This should be set if read by a msh file |
---|
202 | #self.zone = zone |
---|
203 | #self.xllcorner = xllcorner |
---|
204 | #self.yllcorner = yllcorner |
---|
205 | |
---|
206 | |
---|
207 | #Checkpointing and storage |
---|
208 | from anuga.config import default_datadir |
---|
209 | self.datadir = default_datadir |
---|
210 | self.simulation_name = 'domain' |
---|
211 | self.checkpoint = False |
---|
212 | |
---|
213 | #MH310505 To avoid calculating the flux across each edge twice, keep an integer (boolean) array, |
---|
214 | #to be used during the flux calculation |
---|
215 | N = len(self) #number_of_triangles |
---|
216 | self.already_computed_flux = zeros((N, 3), Int) |
---|
217 | |
---|
218 | # Storage for maximal speeds computed for each triangle by compute_fluxes |
---|
219 | # This is used for diagnostics only |
---|
220 | self.max_speed = zeros(N, Float) |
---|
221 | |
---|
222 | if mesh_filename is not None: |
---|
223 | # If the mesh file passed any quantity values |
---|
224 | # , initialise with these values. |
---|
225 | if verbose: print 'Domain: Initialising quantity values' |
---|
226 | self.set_quantity_vertices_dict(vertex_quantity_dict) |
---|
227 | |
---|
228 | |
---|
229 | if verbose: print 'Domain: Done' |
---|
230 | |
---|
231 | |
---|
232 | |
---|
233 | |
---|
234 | def set_default_order(self, n): |
---|
235 | """Set default (spatial) order to either 1 or 2 |
---|
236 | """ |
---|
237 | |
---|
238 | msg = 'Default order must be either 1 or 2. I got %s' %n |
---|
239 | assert n in [1,2], msg |
---|
240 | |
---|
241 | self.default_order = n |
---|
242 | self._order_ = self.default_order |
---|
243 | |
---|
244 | |
---|
245 | #Public interface to Domain |
---|
246 | def get_conserved_quantities(self, vol_id, vertex=None, edge=None): |
---|
247 | """Get conserved quantities at volume vol_id |
---|
248 | |
---|
249 | If vertex is specified use it as index for vertex values |
---|
250 | If edge is specified use it as index for edge values |
---|
251 | If neither are specified use centroid values |
---|
252 | If both are specified an exeception is raised |
---|
253 | |
---|
254 | Return value: Vector of length == number_of_conserved quantities |
---|
255 | |
---|
256 | """ |
---|
257 | |
---|
258 | from Numeric import zeros, Float |
---|
259 | |
---|
260 | if not (vertex is None or edge is None): |
---|
261 | msg = 'Values for both vertex and edge was specified.' |
---|
262 | msg += 'Only one (or none) is allowed.' |
---|
263 | raise msg |
---|
264 | |
---|
265 | q = zeros( len(self.conserved_quantities), Float) |
---|
266 | |
---|
267 | for i, name in enumerate(self.conserved_quantities): |
---|
268 | Q = self.quantities[name] |
---|
269 | if vertex is not None: |
---|
270 | q[i] = Q.vertex_values[vol_id, vertex] |
---|
271 | elif edge is not None: |
---|
272 | q[i] = Q.edge_values[vol_id, edge] |
---|
273 | else: |
---|
274 | q[i] = Q.centroid_values[vol_id] |
---|
275 | |
---|
276 | return q |
---|
277 | |
---|
278 | def set_time(self, time=0.0): |
---|
279 | """Set the model time (seconds)""" |
---|
280 | |
---|
281 | self.time = time |
---|
282 | |
---|
283 | def get_time(self): |
---|
284 | """Get the model time (seconds)""" |
---|
285 | |
---|
286 | return self.time |
---|
287 | |
---|
288 | def set_quantity_vertices_dict(self, quantity_dict): |
---|
289 | """Set values for named quantities. |
---|
290 | The index is the quantity |
---|
291 | |
---|
292 | name: Name of quantity |
---|
293 | X: Compatible list, Numeric array, const or function (see below) |
---|
294 | |
---|
295 | The values will be stored in elements following their |
---|
296 | internal ordering. |
---|
297 | |
---|
298 | """ |
---|
299 | for key in quantity_dict.keys(): |
---|
300 | self.set_quantity(key, quantity_dict[key], location='vertices') |
---|
301 | |
---|
302 | |
---|
303 | def set_quantity(self, name, *args, **kwargs): |
---|
304 | """Set values for named quantity |
---|
305 | |
---|
306 | |
---|
307 | One keyword argument is documented here: |
---|
308 | expression = None, # Arbitrary expression |
---|
309 | |
---|
310 | expression: |
---|
311 | Arbitrary expression involving quantity names |
---|
312 | |
---|
313 | See Quantity.set_values for further documentation. |
---|
314 | """ |
---|
315 | |
---|
316 | #FIXME (Ole): Allow new quantities here |
---|
317 | #from quantity import Quantity, Conserved_quantity |
---|
318 | #Create appropriate quantity object |
---|
319 | ##if name in self.conserved_quantities: |
---|
320 | ## self.quantities[name] = Conserved_quantity(self) |
---|
321 | ##else: |
---|
322 | ## self.quantities[name] = Quantity(self) |
---|
323 | |
---|
324 | |
---|
325 | #Do the expression stuff |
---|
326 | if kwargs.has_key('expression'): |
---|
327 | expression = kwargs['expression'] |
---|
328 | del kwargs['expression'] |
---|
329 | |
---|
330 | Q = self.create_quantity_from_expression(expression) |
---|
331 | kwargs['quantity'] = Q |
---|
332 | |
---|
333 | #Assign values |
---|
334 | self.quantities[name].set_values(*args, **kwargs) |
---|
335 | |
---|
336 | |
---|
337 | def get_quantity_names(self): |
---|
338 | """Get a list of all the quantity names that this domain is aware of. |
---|
339 | Any value in the result should be a valid input to get_quantity. |
---|
340 | """ |
---|
341 | return self.quantities.keys() |
---|
342 | |
---|
343 | def get_quantity(self, name, location='vertices', indices = None): |
---|
344 | """Get quantity object. |
---|
345 | |
---|
346 | name: Name of quantity |
---|
347 | |
---|
348 | See methods inside the quantity object for more options |
---|
349 | |
---|
350 | FIXME: clean input args |
---|
351 | """ |
---|
352 | |
---|
353 | return self.quantities[name] #.get_values( location, indices = indices) |
---|
354 | |
---|
355 | |
---|
356 | def get_quantity_object(self, name): |
---|
357 | """Get object for named quantity |
---|
358 | |
---|
359 | name: Name of quantity |
---|
360 | |
---|
361 | FIXME: Obsolete |
---|
362 | """ |
---|
363 | |
---|
364 | print 'get_quantity_object has been deprecated. Please use get_quantity' |
---|
365 | return self.quantities[name] |
---|
366 | |
---|
367 | |
---|
368 | def create_quantity_from_expression(self, expression): |
---|
369 | """Create new quantity from other quantities using arbitrary expression |
---|
370 | |
---|
371 | Combine existing quantities in domain using expression and return |
---|
372 | result as a new quantity. |
---|
373 | |
---|
374 | Note, the new quantity could e.g. be used in set_quantity |
---|
375 | |
---|
376 | Valid expressions are limited to operators defined in class Quantity |
---|
377 | |
---|
378 | Examples creating derived quantities: |
---|
379 | |
---|
380 | Depth = domain.create_quantity_from_expression('stage-elevation') |
---|
381 | |
---|
382 | exp = '(xmomentum*xmomentum + ymomentum*ymomentum)**0.5') |
---|
383 | Absolute_momentum = domain.create_quantity_from_expression(exp) |
---|
384 | |
---|
385 | """ |
---|
386 | |
---|
387 | from anuga.abstract_2d_finite_volumes.util import\ |
---|
388 | apply_expression_to_dictionary |
---|
389 | |
---|
390 | return apply_expression_to_dictionary(expression, self.quantities) |
---|
391 | |
---|
392 | |
---|
393 | |
---|
394 | #def modify_boundary(self, boundary_map): |
---|
395 | # """Modify existing boundary by elements in boundary map# |
---|
396 | # |
---|
397 | # Input:# |
---|
398 | # |
---|
399 | # boundary_map: Dictionary mapping tags to boundary objects |
---|
400 | # |
---|
401 | # See set_boundary for more details on how this works |
---|
402 | # |
---|
403 | # OBSOLETE |
---|
404 | # """ |
---|
405 | # |
---|
406 | # for key in boundary_map.keys(): |
---|
407 | # self.boundary_map[key] = boundary_map[key] |
---|
408 | # |
---|
409 | # self.set_boundary(self.boundary_map) |
---|
410 | |
---|
411 | |
---|
412 | |
---|
413 | def set_boundary(self, boundary_map): |
---|
414 | """Associate boundary objects with tagged boundary segments. |
---|
415 | |
---|
416 | Input boundary_map is a dictionary of boundary objects keyed |
---|
417 | by symbolic tags to matched against tags in the internal dictionary |
---|
418 | self.boundary. |
---|
419 | |
---|
420 | As result one pointer to a boundary object is stored for each vertex |
---|
421 | in the list self.boundary_objects. |
---|
422 | More entries may point to the same boundary object |
---|
423 | |
---|
424 | Schematically the mapping is from two dictionaries to one list |
---|
425 | where the index is used as pointer to the boundary_values arrays |
---|
426 | within each quantity. |
---|
427 | |
---|
428 | self.boundary: (vol_id, edge_id): tag |
---|
429 | boundary_map (input): tag: boundary_object |
---|
430 | ---------------------------------------------- |
---|
431 | self.boundary_objects: ((vol_id, edge_id), boundary_object) |
---|
432 | |
---|
433 | |
---|
434 | Pre-condition: |
---|
435 | self.boundary has been built. |
---|
436 | |
---|
437 | Post-condition: |
---|
438 | self.boundary_objects is built |
---|
439 | |
---|
440 | If a tag from the domain doesn't appear in the input dictionary an |
---|
441 | exception is raised. |
---|
442 | However, if a tag is not used to the domain, no error is thrown. |
---|
443 | FIXME: This would lead to implementation of a |
---|
444 | default boundary condition |
---|
445 | |
---|
446 | Note: If a segment is listed in the boundary dictionary and if it is |
---|
447 | not None, it *will* become a boundary - |
---|
448 | even if there is a neighbouring triangle. |
---|
449 | This would be the case for internal boundaries |
---|
450 | |
---|
451 | Boundary objects that are None will be skipped. |
---|
452 | |
---|
453 | If a boundary_map has already been set |
---|
454 | (i.e. set_boundary has been called before), the old boundary map |
---|
455 | will be updated with new values. The new map need not define all |
---|
456 | boundary tags, and can thus change only those that are needed. |
---|
457 | |
---|
458 | FIXME: If set_boundary is called multiple times and if Boundary |
---|
459 | object is changed into None, the neighbour structure will not be |
---|
460 | restored!!! |
---|
461 | |
---|
462 | |
---|
463 | """ |
---|
464 | |
---|
465 | if self.boundary_map is None: |
---|
466 | # This the first call to set_boundary. Store |
---|
467 | # map for later updates and for use with boundary_stats. |
---|
468 | self.boundary_map = boundary_map |
---|
469 | else: |
---|
470 | # This is a modification of an already existing map |
---|
471 | # Update map an proceed normally |
---|
472 | |
---|
473 | for key in boundary_map.keys(): |
---|
474 | self.boundary_map[key] = boundary_map[key] |
---|
475 | |
---|
476 | |
---|
477 | #FIXME: Try to remove the sorting and fix test_mesh.py |
---|
478 | x = self.boundary.keys() |
---|
479 | x.sort() |
---|
480 | |
---|
481 | #Loop through edges that lie on the boundary and associate them with |
---|
482 | #callable boundary objects depending on their tags |
---|
483 | self.boundary_objects = [] |
---|
484 | for k, (vol_id, edge_id) in enumerate(x): |
---|
485 | tag = self.boundary[ (vol_id, edge_id) ] |
---|
486 | |
---|
487 | if self.boundary_map.has_key(tag): |
---|
488 | B = self.boundary_map[tag] #Get callable boundary object |
---|
489 | |
---|
490 | if B is not None: |
---|
491 | self.boundary_objects.append( ((vol_id, edge_id), B) ) |
---|
492 | self.neighbours[vol_id, edge_id] = \ |
---|
493 | -len(self.boundary_objects) |
---|
494 | else: |
---|
495 | pass |
---|
496 | #FIXME: Check and perhaps fix neighbour structure |
---|
497 | |
---|
498 | |
---|
499 | else: |
---|
500 | msg = 'ERROR (domain.py): Tag "%s" has not been ' %tag |
---|
501 | msg += 'bound to a boundary object.\n' |
---|
502 | msg += 'All boundary tags defined in domain must appear ' |
---|
503 | msg += 'in the supplied dictionary.\n' |
---|
504 | msg += 'The tags are: %s' %self.get_boundary_tags() |
---|
505 | raise msg |
---|
506 | |
---|
507 | |
---|
508 | def set_region(self, *args, **kwargs): |
---|
509 | """ |
---|
510 | This method is used to set quantities based on a regional tag. |
---|
511 | |
---|
512 | It is most often called with the following parameters; |
---|
513 | (self, tag, quantity, X, location='vertices') |
---|
514 | tag: the name of the regional tag used to specify the region |
---|
515 | quantity: Name of quantity to change |
---|
516 | X: const or function - how the quantity is changed |
---|
517 | location: Where values are to be stored. |
---|
518 | Permissible options are: vertices, centroid and unique vertices |
---|
519 | |
---|
520 | A callable region class or a list of callable region classes |
---|
521 | can also be passed into this function. |
---|
522 | """ |
---|
523 | #print "*args", args |
---|
524 | #print "**kwargs", kwargs |
---|
525 | if len(args) == 1: |
---|
526 | self._set_region(*args, **kwargs) |
---|
527 | else: |
---|
528 | #Assume it is arguments for the region.set_region function |
---|
529 | func = region_set_region(*args, **kwargs) |
---|
530 | self._set_region(func) |
---|
531 | |
---|
532 | |
---|
533 | def _set_region(self, functions): |
---|
534 | # The order of functions in the list is used. |
---|
535 | if type(functions) not in [types.ListType,types.TupleType]: |
---|
536 | functions = [functions] |
---|
537 | for function in functions: |
---|
538 | for tag in self.tagged_elements.keys(): |
---|
539 | function(tag, self.tagged_elements[tag], self) |
---|
540 | |
---|
541 | |
---|
542 | #MISC |
---|
543 | def check_integrity(self): |
---|
544 | Mesh.check_integrity(self) |
---|
545 | |
---|
546 | for quantity in self.conserved_quantities: |
---|
547 | msg = 'Conserved quantities must be a subset of all quantities' |
---|
548 | assert quantity in self.quantities, msg |
---|
549 | |
---|
550 | ##assert hasattr(self, 'boundary_objects') |
---|
551 | |
---|
552 | def write_time(self, track_speeds=False): |
---|
553 | print self.timestepping_statistics(track_speeds) |
---|
554 | |
---|
555 | |
---|
556 | def timestepping_statistics(self, track_speeds=False): |
---|
557 | """Return string with time stepping statistics for printing or logging |
---|
558 | |
---|
559 | Optional boolean keyword track_speeds decides whether to report location of |
---|
560 | smallest timestep as well as a histogram and percentile report. |
---|
561 | """ |
---|
562 | |
---|
563 | from anuga.utilities.numerical_tools import histogram, create_bins |
---|
564 | |
---|
565 | |
---|
566 | msg = '' |
---|
567 | if self.min_timestep == self.max_timestep: |
---|
568 | msg += 'Time = %.4f, delta t = %.8f, steps=%d (%d)'\ |
---|
569 | %(self.time, self.min_timestep, self.number_of_steps, |
---|
570 | self.number_of_first_order_steps) |
---|
571 | elif self.min_timestep > self.max_timestep: |
---|
572 | msg += 'Time = %.4f, steps=%d (%d)'\ |
---|
573 | %(self.time, self.number_of_steps, |
---|
574 | self.number_of_first_order_steps) |
---|
575 | else: |
---|
576 | msg += 'Time = %.4f, delta t in [%.8f, %.8f], steps=%d (%d)'\ |
---|
577 | %(self.time, self.min_timestep, |
---|
578 | self.max_timestep, self.number_of_steps, |
---|
579 | self.number_of_first_order_steps) |
---|
580 | |
---|
581 | if track_speeds is True: |
---|
582 | msg += '\n' |
---|
583 | |
---|
584 | |
---|
585 | #Setup 10 bins for speed histogram |
---|
586 | bins = create_bins(self.max_speed, 10) |
---|
587 | hist = histogram(self.max_speed, bins) |
---|
588 | |
---|
589 | msg += '------------------------------------------------\n' |
---|
590 | msg += ' Speeds in [%f, %f]\n' %(min(self.max_speed), max(self.max_speed)) |
---|
591 | msg += ' Histogram:\n' |
---|
592 | |
---|
593 | hi = bins[0] |
---|
594 | for i, count in enumerate(hist): |
---|
595 | lo = hi |
---|
596 | if i+1 < len(bins): |
---|
597 | #Open upper interval |
---|
598 | hi = bins[i+1] |
---|
599 | msg += ' [%f, %f[: %d\n' %(lo, hi, count) |
---|
600 | else: |
---|
601 | #Closed upper interval |
---|
602 | hi = max(self.max_speed) |
---|
603 | msg += ' [%f, %f]: %d\n' %(lo, hi, count) |
---|
604 | |
---|
605 | |
---|
606 | N = len(self.max_speed) |
---|
607 | if N > 10: |
---|
608 | msg += ' Percentiles (10%):\n' |
---|
609 | speed = self.max_speed.tolist() |
---|
610 | speed.sort() |
---|
611 | |
---|
612 | k = 0 |
---|
613 | lower = min(speed) |
---|
614 | for i, a in enumerate(speed): |
---|
615 | if i % (N/10) == 0 and i != 0: #For every 10% of the sorted speeds |
---|
616 | msg += ' %d speeds in [%f, %f]\n' %(i-k, lower, a) |
---|
617 | lower = a |
---|
618 | k = i |
---|
619 | |
---|
620 | msg += ' %d speeds in [%f, %f]\n'\ |
---|
621 | %(N-k, lower, max(speed)) |
---|
622 | |
---|
623 | |
---|
624 | |
---|
625 | |
---|
626 | |
---|
627 | |
---|
628 | # Find index of largest computed flux speed |
---|
629 | k = argmax(self.max_speed) |
---|
630 | |
---|
631 | x, y = self.get_centroid_coordinates()[k] |
---|
632 | |
---|
633 | msg += ' Triangle #%d with centroid (%.4f, %.4f) ' %(k, x, y) |
---|
634 | msg += 'had the largest computed speed: %.4f m/s\n' %(self.max_speed[k]) |
---|
635 | |
---|
636 | # Report all quantity values at vertices |
---|
637 | msg += ' Quantity \t vertex values\t\t\t\t\t centroid values\n' |
---|
638 | for name in self.quantities: |
---|
639 | q = self.quantities[name] |
---|
640 | X,Y,A,V = q.get_vertex_values() |
---|
641 | |
---|
642 | s = ' %s:\t %.4f,\t %.4f,\t %.4f,\t %.4f\n'\ |
---|
643 | %(name, A[3*k], A[3*k+1], A[3*k+2], q.get_values(location='centroids')[k]) |
---|
644 | |
---|
645 | msg += s |
---|
646 | |
---|
647 | return msg |
---|
648 | |
---|
649 | |
---|
650 | def write_boundary_statistics(self, quantities = None, tags = None): |
---|
651 | print self.boundary_statistics(quantities, tags) |
---|
652 | |
---|
653 | def boundary_statistics(self, quantities = None, tags = None): |
---|
654 | """Output statistics about boundary forcing at each timestep |
---|
655 | |
---|
656 | |
---|
657 | Input: |
---|
658 | quantities: either None, a string or a list of strings naming the quantities to be reported |
---|
659 | tags: either None, a string or a list of strings naming the tags to be reported |
---|
660 | |
---|
661 | |
---|
662 | Example output: |
---|
663 | Tag 'wall': |
---|
664 | stage in [2, 5.5] |
---|
665 | xmomentum in [] |
---|
666 | ymomentum in [] |
---|
667 | Tag 'ocean' |
---|
668 | |
---|
669 | |
---|
670 | If quantities are specified only report on those. Otherwise take all conserved quantities. |
---|
671 | If tags are specified only report on those, otherwise take all tags. |
---|
672 | |
---|
673 | """ |
---|
674 | |
---|
675 | #Input checks |
---|
676 | import types, string |
---|
677 | |
---|
678 | if quantities is None: |
---|
679 | quantities = self.conserved_quantities |
---|
680 | elif type(quantities) == types.StringType: |
---|
681 | quantities = [quantities] #Turn it into a list |
---|
682 | |
---|
683 | msg = 'Keyword argument quantities must be either None, ' |
---|
684 | msg += 'string or list. I got %s' %str(quantities) |
---|
685 | assert type(quantities) == types.ListType, msg |
---|
686 | |
---|
687 | |
---|
688 | if tags is None: |
---|
689 | tags = self.get_boundary_tags() |
---|
690 | elif type(tags) == types.StringType: |
---|
691 | tags = [tags] #Turn it into a list |
---|
692 | |
---|
693 | msg = 'Keyword argument tags must be either None, ' |
---|
694 | msg += 'string or list. I got %s' %str(tags) |
---|
695 | assert type(tags) == types.ListType, msg |
---|
696 | |
---|
697 | #Determine width of longest quantity name (for cosmetic purposes) |
---|
698 | maxwidth = 0 |
---|
699 | for name in quantities: |
---|
700 | w = len(name) |
---|
701 | if w > maxwidth: |
---|
702 | maxwidth = w |
---|
703 | |
---|
704 | #Output stats |
---|
705 | msg = 'Boundary values at time %.4f:\n' %self.time |
---|
706 | for tag in tags: |
---|
707 | msg += ' %s:\n' %tag |
---|
708 | |
---|
709 | for name in quantities: |
---|
710 | q = self.quantities[name] |
---|
711 | |
---|
712 | #Find range of boundary values for tag and q |
---|
713 | maxval = minval = None |
---|
714 | for i, ((vol_id, edge_id), B) in\ |
---|
715 | enumerate(self.boundary_objects): |
---|
716 | if self.boundary[(vol_id, edge_id)] == tag: |
---|
717 | v = q.boundary_values[i] |
---|
718 | if minval is None or v < minval: minval = v |
---|
719 | if maxval is None or v > maxval: maxval = v |
---|
720 | |
---|
721 | if minval is None or maxval is None: |
---|
722 | msg += ' Sorry no information available about' +\ |
---|
723 | ' tag %s and quantity %s\n' %(tag, name) |
---|
724 | else: |
---|
725 | msg += ' %s in [%12.8f, %12.8f]\n'\ |
---|
726 | %(string.ljust(name, maxwidth), minval, maxval) |
---|
727 | |
---|
728 | |
---|
729 | return msg |
---|
730 | |
---|
731 | |
---|
732 | def get_name(self): |
---|
733 | return self.simulation_name |
---|
734 | |
---|
735 | def set_name(self, name): |
---|
736 | """Assign a name to this simulation. |
---|
737 | This will be used to identify the output sww file. |
---|
738 | |
---|
739 | """ |
---|
740 | if name.endswith('.sww'): |
---|
741 | name = name[:-4] |
---|
742 | |
---|
743 | self.simulation_name = name |
---|
744 | |
---|
745 | def get_datadir(self): |
---|
746 | return self.datadir |
---|
747 | |
---|
748 | def set_datadir(self, name): |
---|
749 | self.datadir = name |
---|
750 | |
---|
751 | |
---|
752 | |
---|
753 | #def set_defaults(self): |
---|
754 | # """Set default values for uninitialised quantities. |
---|
755 | # Should be overridden or specialised by specific modules |
---|
756 | # """# |
---|
757 | # |
---|
758 | # for name in self.conserved_quantities + self.other_quantities: |
---|
759 | # self.set_quantity(name, 0.0) |
---|
760 | |
---|
761 | |
---|
762 | ########################### |
---|
763 | #Main components of evolve |
---|
764 | |
---|
765 | def evolve(self, |
---|
766 | yieldstep = None, |
---|
767 | finaltime = None, |
---|
768 | duration = None, |
---|
769 | skip_initial_step = False): |
---|
770 | """Evolve model through time starting from self.starttime. |
---|
771 | |
---|
772 | |
---|
773 | yieldstep: Interval between yields where results are stored, |
---|
774 | statistics written and domain inspected or |
---|
775 | possibly modified. If omitted the internal predefined |
---|
776 | max timestep is used. |
---|
777 | Internally, smaller timesteps may be taken. |
---|
778 | |
---|
779 | duration: Duration of simulation |
---|
780 | |
---|
781 | finaltime: Time where simulation should end |
---|
782 | |
---|
783 | If both duration and finaltime are given an exception is thrown. |
---|
784 | |
---|
785 | |
---|
786 | skip_initial_step: Boolean flag that decides whether the first |
---|
787 | yield step is skipped or not. This is useful for example to avoid |
---|
788 | duplicate steps when multiple evolve processes are dove tailed. |
---|
789 | |
---|
790 | |
---|
791 | Evolve is implemented as a generator and is to be called as such, e.g. |
---|
792 | |
---|
793 | for t in domain.evolve(yieldstep, finaltime): |
---|
794 | <Do something with domain and t> |
---|
795 | |
---|
796 | |
---|
797 | All times are given in seconds |
---|
798 | |
---|
799 | """ |
---|
800 | |
---|
801 | from anuga.config import min_timestep, max_timestep, epsilon |
---|
802 | |
---|
803 | #FIXME: Maybe lump into a larger check prior to evolving |
---|
804 | msg = 'Boundary tags must be bound to boundary objects before ' |
---|
805 | msg += 'evolving system, ' |
---|
806 | msg += 'e.g. using the method set_boundary.\n' |
---|
807 | msg += 'This system has the boundary tags %s '\ |
---|
808 | %self.get_boundary_tags() |
---|
809 | assert hasattr(self, 'boundary_objects'), msg |
---|
810 | |
---|
811 | |
---|
812 | if yieldstep is None: |
---|
813 | yieldstep = max_timestep |
---|
814 | else: |
---|
815 | yieldstep = float(yieldstep) |
---|
816 | |
---|
817 | self._order_ = self.default_order |
---|
818 | |
---|
819 | |
---|
820 | if finaltime is not None and duration is not None: |
---|
821 | print 'F', finaltime, duration |
---|
822 | msg = 'Only one of finaltime and duration may be specified' |
---|
823 | raise msg |
---|
824 | else: |
---|
825 | if finaltime is not None: |
---|
826 | self.finaltime = float(finaltime) |
---|
827 | if duration is not None: |
---|
828 | self.finaltime = self.starttime + float(duration) |
---|
829 | |
---|
830 | |
---|
831 | |
---|
832 | |
---|
833 | self.yieldtime = 0.0 #Time between 'yields' |
---|
834 | |
---|
835 | #Initialise interval of timestep sizes (for reporting only) |
---|
836 | self.min_timestep = max_timestep |
---|
837 | self.max_timestep = min_timestep |
---|
838 | self.number_of_steps = 0 |
---|
839 | self.number_of_first_order_steps = 0 |
---|
840 | |
---|
841 | #update ghosts |
---|
842 | self.update_ghosts() |
---|
843 | |
---|
844 | #Initial update of vertex and edge values |
---|
845 | self.distribute_to_vertices_and_edges() |
---|
846 | |
---|
847 | #Initial update boundary values |
---|
848 | self.update_boundary() |
---|
849 | |
---|
850 | #Or maybe restore from latest checkpoint |
---|
851 | if self.checkpoint is True: |
---|
852 | self.goto_latest_checkpoint() |
---|
853 | |
---|
854 | if skip_initial_step is False: |
---|
855 | yield(self.time) #Yield initial values |
---|
856 | |
---|
857 | while True: |
---|
858 | |
---|
859 | #Compute fluxes across each element edge |
---|
860 | self.compute_fluxes() |
---|
861 | |
---|
862 | #Update timestep to fit yieldstep and finaltime |
---|
863 | self.update_timestep(yieldstep, finaltime) |
---|
864 | |
---|
865 | #Update conserved quantities |
---|
866 | self.update_conserved_quantities() |
---|
867 | |
---|
868 | #update ghosts |
---|
869 | self.update_ghosts() |
---|
870 | |
---|
871 | #Update vertex and edge values |
---|
872 | self.distribute_to_vertices_and_edges() |
---|
873 | |
---|
874 | #Update boundary values |
---|
875 | self.update_boundary() |
---|
876 | |
---|
877 | #Update time |
---|
878 | self.time += self.timestep |
---|
879 | self.yieldtime += self.timestep |
---|
880 | self.number_of_steps += 1 |
---|
881 | if self._order_ == 1: |
---|
882 | self.number_of_first_order_steps += 1 |
---|
883 | |
---|
884 | #Yield results |
---|
885 | if finaltime is not None and self.time >= finaltime-epsilon: |
---|
886 | |
---|
887 | if self.time > finaltime: |
---|
888 | #FIXME (Ole, 30 April 2006): Do we need this check? |
---|
889 | print 'WARNING (domain.py): time overshot finaltime. Contact Ole.Nielsen@ga.gov.au' |
---|
890 | self.time = finaltime |
---|
891 | |
---|
892 | # Yield final time and stop |
---|
893 | self.time = finaltime |
---|
894 | yield(self.time) |
---|
895 | break |
---|
896 | |
---|
897 | |
---|
898 | if self.yieldtime >= yieldstep: |
---|
899 | # Yield (intermediate) time and allow inspection of domain |
---|
900 | |
---|
901 | if self.checkpoint is True: |
---|
902 | self.store_checkpoint() |
---|
903 | self.delete_old_checkpoints() |
---|
904 | |
---|
905 | # Pass control on to outer loop for more specific actions |
---|
906 | yield(self.time) |
---|
907 | |
---|
908 | # Reinitialise |
---|
909 | self.yieldtime = 0.0 |
---|
910 | self.min_timestep = max_timestep |
---|
911 | self.max_timestep = min_timestep |
---|
912 | self.number_of_steps = 0 |
---|
913 | self.number_of_first_order_steps = 0 |
---|
914 | |
---|
915 | |
---|
916 | def evolve_to_end(self, finaltime = 1.0): |
---|
917 | """Iterate evolve all the way to the end |
---|
918 | """ |
---|
919 | |
---|
920 | for _ in self.evolve(yieldstep=None, finaltime=finaltime): |
---|
921 | pass |
---|
922 | |
---|
923 | |
---|
924 | |
---|
925 | def update_boundary(self): |
---|
926 | """Go through list of boundary objects and update boundary values |
---|
927 | for all conserved quantities on boundary. |
---|
928 | """ |
---|
929 | |
---|
930 | #FIXME: Update only those that change (if that can be worked out) |
---|
931 | #FIXME: Boundary objects should not include ghost nodes. |
---|
932 | for i, ((vol_id, edge_id), B) in enumerate(self.boundary_objects): |
---|
933 | if B is None: |
---|
934 | print 'WARNING: Ignored boundary segment %d (None)' |
---|
935 | else: |
---|
936 | q = B.evaluate(vol_id, edge_id) |
---|
937 | |
---|
938 | for j, name in enumerate(self.conserved_quantities): |
---|
939 | Q = self.quantities[name] |
---|
940 | Q.boundary_values[i] = q[j] |
---|
941 | |
---|
942 | |
---|
943 | def compute_fluxes(self): |
---|
944 | msg = 'Method compute_fluxes must be overridden by Domain subclass' |
---|
945 | raise msg |
---|
946 | |
---|
947 | |
---|
948 | def update_timestep(self, yieldstep, finaltime): |
---|
949 | |
---|
950 | from anuga.config import min_timestep, max_timestep |
---|
951 | |
---|
952 | # self.timestep is calculated from speed of characteristics |
---|
953 | # Apply CFL condition here |
---|
954 | timestep = min(self.CFL*self.timestep,max_timestep) |
---|
955 | |
---|
956 | #Record maximal and minimal values of timestep for reporting |
---|
957 | self.max_timestep = max(timestep, self.max_timestep) |
---|
958 | self.min_timestep = min(timestep, self.min_timestep) |
---|
959 | |
---|
960 | #Protect against degenerate time steps |
---|
961 | if timestep < min_timestep: |
---|
962 | |
---|
963 | #Number of consecutive small steps taken b4 taking action |
---|
964 | self.smallsteps += 1 |
---|
965 | |
---|
966 | if self.smallsteps > self.max_smallsteps: |
---|
967 | self.smallsteps = 0 #Reset |
---|
968 | |
---|
969 | if self._order_ == 1: |
---|
970 | msg = 'WARNING: Too small timestep %.16f reached '\ |
---|
971 | %timestep |
---|
972 | msg += 'even after %d steps of 1 order scheme'\ |
---|
973 | %self.max_smallsteps |
---|
974 | print msg |
---|
975 | timestep = min_timestep #Try enforcing min_step |
---|
976 | |
---|
977 | #raise msg |
---|
978 | else: |
---|
979 | #Try to overcome situation by switching to 1 order |
---|
980 | self._order_ = 1 |
---|
981 | |
---|
982 | else: |
---|
983 | self.smallsteps = 0 |
---|
984 | if self._order_ == 1 and self.default_order == 2: |
---|
985 | self._order_ = 2 |
---|
986 | |
---|
987 | |
---|
988 | #Ensure that final time is not exceeded |
---|
989 | if finaltime is not None and self.time + timestep > finaltime : |
---|
990 | timestep = finaltime-self.time |
---|
991 | |
---|
992 | #Ensure that model time is aligned with yieldsteps |
---|
993 | if self.yieldtime + timestep > yieldstep: |
---|
994 | timestep = yieldstep-self.yieldtime |
---|
995 | |
---|
996 | self.timestep = timestep |
---|
997 | |
---|
998 | |
---|
999 | |
---|
1000 | def compute_forcing_terms(self): |
---|
1001 | """If there are any forcing functions driving the system |
---|
1002 | they should be defined in Domain subclass and appended to |
---|
1003 | the list self.forcing_terms |
---|
1004 | """ |
---|
1005 | |
---|
1006 | for f in self.forcing_terms: |
---|
1007 | f(self) |
---|
1008 | |
---|
1009 | |
---|
1010 | |
---|
1011 | def update_conserved_quantities(self): |
---|
1012 | """Update vectors of conserved quantities using previously |
---|
1013 | computed fluxes specified forcing functions. |
---|
1014 | """ |
---|
1015 | |
---|
1016 | from Numeric import ones, sum, equal, Float |
---|
1017 | |
---|
1018 | N = len(self) #number_of_triangles |
---|
1019 | d = len(self.conserved_quantities) |
---|
1020 | |
---|
1021 | timestep = self.timestep |
---|
1022 | |
---|
1023 | #Compute forcing terms |
---|
1024 | self.compute_forcing_terms() |
---|
1025 | |
---|
1026 | #Update conserved_quantities |
---|
1027 | for name in self.conserved_quantities: |
---|
1028 | Q = self.quantities[name] |
---|
1029 | Q.update(timestep) |
---|
1030 | |
---|
1031 | #Clean up |
---|
1032 | #Note that Q.explicit_update is reset by compute_fluxes |
---|
1033 | |
---|
1034 | #MH090605 commented out the following since semi_implicit_update is now re-initialized |
---|
1035 | #at the end of the _update function in quantity_ext.c (This is called by the |
---|
1036 | #preceeding Q.update(timestep) statement above). |
---|
1037 | #For run_profile.py with N=128, the time of update_conserved_quantities is cut from 14.00 secs |
---|
1038 | #to 8.35 secs |
---|
1039 | |
---|
1040 | #Q.semi_implicit_update[:] = 0.0 |
---|
1041 | |
---|
1042 | def update_ghosts(self): |
---|
1043 | pass |
---|
1044 | |
---|
1045 | def distribute_to_vertices_and_edges(self): |
---|
1046 | """Extrapolate conserved quantities from centroid to |
---|
1047 | vertices and edge-midpoints for each volume |
---|
1048 | |
---|
1049 | Default implementation is straight first order, |
---|
1050 | i.e. constant values throughout each element and |
---|
1051 | no reference to non-conserved quantities. |
---|
1052 | """ |
---|
1053 | |
---|
1054 | for name in self.conserved_quantities: |
---|
1055 | Q = self.quantities[name] |
---|
1056 | if self._order_ == 1: |
---|
1057 | Q.extrapolate_first_order() |
---|
1058 | elif self._order_ == 2: |
---|
1059 | Q.extrapolate_second_order() |
---|
1060 | Q.limit() |
---|
1061 | else: |
---|
1062 | raise 'Unknown order' |
---|
1063 | Q.interpolate_from_vertices_to_edges() |
---|
1064 | |
---|
1065 | |
---|
1066 | def centroid_norm(self, quantity, normfunc): |
---|
1067 | """Calculate the norm of the centroid values |
---|
1068 | of a specific quantity, using normfunc. |
---|
1069 | |
---|
1070 | normfunc should take a list to a float. |
---|
1071 | |
---|
1072 | common normfuncs are provided in the module utilities.norms |
---|
1073 | """ |
---|
1074 | return normfunc(self.quantities[quantity].centroid_values) |
---|
1075 | |
---|
1076 | |
---|
1077 | |
---|
1078 | ############################################## |
---|
1079 | #Initialise module |
---|
1080 | |
---|
1081 | #Optimisation with psyco |
---|
1082 | from anuga.config import use_psyco |
---|
1083 | if use_psyco: |
---|
1084 | try: |
---|
1085 | import psyco |
---|
1086 | except: |
---|
1087 | import os |
---|
1088 | if os.name == 'posix' and os.uname()[4] == 'x86_64': |
---|
1089 | pass |
---|
1090 | #Psyco isn't supported on 64 bit systems, but it doesn't matter |
---|
1091 | else: |
---|
1092 | msg = 'WARNING: psyco (speedup) could not import'+\ |
---|
1093 | ', you may want to consider installing it' |
---|
1094 | print msg |
---|
1095 | else: |
---|
1096 | psyco.bind(Domain.update_boundary) |
---|
1097 | #psyco.bind(Domain.update_timestep) #Not worth it |
---|
1098 | psyco.bind(Domain.update_conserved_quantities) |
---|
1099 | psyco.bind(Domain.distribute_to_vertices_and_edges) |
---|
1100 | |
---|
1101 | |
---|
1102 | if __name__ == "__main__": |
---|
1103 | pass |
---|