hedge.discretization

1 # -*- coding: utf8 -*- 2 3 """Global function space discretization.""" 4 5 from __future__ import division 6 7 __copyright__ = "Copyright (C) 2007 Andreas Kloeckner" 8 9 __license__ = """ 10 This program is free software: you can redistribute it and/or modify 11 it under the terms of the GNU General Public License as published by 12 the Free Software Foundation, either version 3 of the License, or 13 (at your option) any later version. 14 15 This program is distributed in the hope that it will be useful, 16 but WITHOUT ANY WARRANTY; without even the implied warranty of 17 MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the 18 GNU General Public License for more details. 19 20 You should have received a copy of the GNU General Public License 21 along with this program. If not, see U{http://www.gnu.org/licenses/}. 22 """ 23 24 25 26 27 import numpy 28 import numpy.linalg as la 29 import hedge.tools 30 import hedge.mesh 31 import hedge.optemplate 32 import hedge._internal 33 from pytools import memoize_method

34 35 36 37 38 -class _FaceGroup(hedge._internal.FaceGroup):

39 - def __init__(self, double_sided, debug):

40 hedge._internal.FaceGroup.__init__(self, double_sided) 41 from hedge.tools import IndexListRegistry 42 self.fil_registry = IndexListRegistry(debug)

43

44 - def register_face_index_list(self, identifier, generator):

45 return self.fil_registry.register(identifier, generator)

46

47 - def commit(self, discr, ldis_loc, ldis_opp):

48 if self.fil_registry.index_lists: 49 self.index_lists = numpy.array( 50 self.fil_registry.index_lists, 51 dtype=numpy.uint32, order="C") 52 del self.fil_registry 53 54 if ldis_loc is None: 55 self.face_count = 0 56 else: 57 self.face_count = ldis_loc.face_count() 58 59 # number elements locally 60 used_bases_and_els = list(set( 61 (side.el_base_index, side.element_id) 62 for fp in self.face_pairs 63 for side in [fp.loc, fp.opp] 64 if side.element_id != hedge._internal.INVALID_ELEMENT)) 65 66 used_bases_and_els.sort() 67 el_id_to_local_number = dict( 68 (bae[1], i) for i, bae in enumerate(used_bases_and_els)) 69 self.local_el_to_global_el_base = numpy.fromiter( 70 (bae[0] for bae in used_bases_and_els), dtype=numpy.uint32) 71 72 for fp in self.face_pairs: 73 for side in [fp.loc, fp.opp]: 74 if side.element_id != hedge._internal.INVALID_ELEMENT: 75 side.local_el_number = el_id_to_local_number[side.element_id] 76 77 # transfer inverse jacobians 78 self.local_el_inverse_jacobians = numpy.fromiter( 79 (abs(discr.mesh.elements[bae[1]].inverse_map.jacobian()) 80 for bae in used_bases_and_els), 81 dtype=float) 82 83 self.ldis_loc = ldis_loc 84 self.ldis_opp = ldis_opp

85

86 87 88 89 -class _ElementGroup(object):

90 """Once fully filled, this structure has the following data members: 91 92 @ivar members: a list of L{hedge.mesh.Element} instances in this group. 93 @ivar member_nrs: a list of the element ID numbers in this group. 94 @ivar local_discretization: an instance of L{hedge.element.Element}. 95 @ivar ranges: a list of C{slice} objects indicating the DOF numbers for 96 each element. Note: This is actually a C++ ElementRanges object. 97 @ivar mass_matrix: The element-local mass matrix M{M}. 98 @ivar inverse_mass_matrix: the element-local inverese mass matrix M{M^{-1}}. 99 @ivar differentiation_matrices: local differentiation matrices M{D_r, D_s, D_t}, 100 i.e. differentiation by M{r, s, t, ....}. 101 @ivar stiffness_matrices: the element-local stiffness matrices M{M*D_r, M*D_s,...}. 102 @ivar jacobians: list of jacobians over all elements 103 @ivar inverse_jacobians: inverses of L{jacobians}. 104 @ivar diff_coefficients: a M{(d,d)}-matrix of coefficient vectors to turn 105 M{(r,s,t)}-differentiation into M{(x,y,z)}. 106 """ 107 pass

108

109 110 111 112 -class _Boundary(object):

113 - def __init__(self, nodes, ranges, vol_indices, face_groups, 114 el_face_to_face_group_and_face_pair={}):

115 self.nodes = nodes 116 self.ranges = ranges 117 self.vol_indices = vol_indices 118 self.face_groups = face_groups 119 self.el_face_to_face_group_and_face_pair = \ 120 el_face_to_face_group_and_face_pair

121

122 - def find_facepair(self, el_face):

123 fg, fp_idx = self.el_face_to_face_group_and_face_pair[el_face] 124 125 return fg.face_pairs[fp_idx]

126

127 - def find_facepair_side(self, el_face):

128 fp = self.find_facepair(el_face) 129 el, face_nbr = el_face 130 131 for flux_face in [fp.loc, fp.opp]: 132 if flux_face.element_id == el.id and flux_face.face_id == face_nbr: 133 return flux_face 134 raise KeyError, "flux face not found in boundary"

135

136 137 138 139 -class OpTemplateFunction:

140 - def __init__(self, discr, pp_optemplate):

141 self.discr = discr 142 self.pp_optemplate = pp_optemplate

143

144 - def __call__(self, **vars):

145 return self.discr.run_preprocessed_optemplate(self.pp_optemplate, vars)

146

147 148 149 150 -class _PointEvaluator(object):

151 - def __init__(self, discr, el_range, interp_coeff):

152 self.discr = discr 153 self.el_range = el_range 154 self.interp_coeff = interp_coeff

155

156 - def __call__(self, field):

157 from hedge.tools import log_shape 158 ls = log_shape(field) 159 if ls != (): 160 result = numpy.zeros(ls, dtype=self.discr.default_scalar_type) 161 from pytools import indices_in_shape 162 for i in indices_in_shape(ls): 163 result[i] = numpy.dot( 164 self.interp_coeff, field[i][self.el_range]) 165 return result 166 else: 167 return numpy.dot(self.interp_coeff, field[self.el_range])

168

169 170 -class Discretization(object):

171 """The global approximation space. 172 173 Instances of this class tie together a local discretization (i.e. polynomials 174 on an elemnent) into a function space on a mesh. They provide creation 175 functions such as interpolating given functions, differential operators and 176 flux lifting operators. 177 """ 178 179 @classmethod

180 - def all_debug_flags(cls):

181 return set([ 182 "ilist_generation", 183 "node_permutation", 184 "print_op_code", 185 "dump_dataflow_graph", 186 ])

187 188 @classmethod

189 - def noninteractive_debug_flags(cls):

190 return set([ 191 "ilist_generation", 192 "node_permutation", 193 ])

194 195 @staticmethod

196 - def get_local_discretization(mesh, local_discretization=None, order=None):

197 if local_discretization is None and order is None: 198 raise ValueError, "must supply either local_discretization or order" 199 if local_discretization is not None and order is not None: 200 raise ValueError, "must supply only one of local_discretization and order" 201 if local_discretization is None: 202 from hedge.element import ELEMENTS 203 from pytools import one 204 ldis_class = one( 205 ldis_class for ldis_class in ELEMENTS 206 if isinstance(mesh.elements[0], ldis_class.geometry)) 207 return ldis_class(order) 208 else: 209 return local_discretization

210

211 - def __init__(self, mesh, local_discretization=None, 212 order=None, debug=set(), default_scalar_type=numpy.float64, 213 run_context=None):

214 """ 215 216 @arg debug: A set of strings indicating which debug checks should 217 be activated. See validity check below for the currently defined 218 set of debug flags. 219 """ 220 221 self.mesh = mesh 222 223 local_discretization = self.get_local_discretization( 224 mesh, local_discretization, order) 225 226 self.dimensions = local_discretization.dimensions 227 228 debug = set(debug) 229 assert not debug.difference(self.all_debug_flags()), "Invalid debug flag specified" 230 self.debug = debug 231 232 self._build_element_groups_and_nodes(local_discretization) 233 self._calculate_local_matrices() 234 self._build_interior_face_groups() 235 236 self.instrumented = False 237 238 self.default_scalar_type = default_scalar_type 239 240 self.exec_functions = {}

241

242 - def close(self):

243 pass

244 245 # instrumentation ---------------------------------------------------------

246 - def create_op_timers(self):

247 from pytools.log import IntervalTimer 248 249 self.gather_timer = IntervalTimer("t_gather", 250 "Time spent gathering fluxes") 251 self.lift_timer = IntervalTimer("t_lift", 252 "Time spent lifting fluxes") 253 self.mass_timer = IntervalTimer("t_mass", 254 "Time spent applying mass operators") 255 self.diff_timer = IntervalTimer("t_diff", 256 "Time spent applying applying differentiation operators") 257 self.vector_math_timer = IntervalTimer("t_vector_math", 258 "Time spent doing vector math") 259 260 return [self.gather_timer, 261 self.lift_timer, 262 self.mass_timer, 263 self.diff_timer, 264 self.vector_math_timer]

265

266 - def add_instrumentation(self, mgr):

267 from pytools.log import IntervalTimer, EventCounter 268 269 self.gather_counter = EventCounter("n_gather", 270 "Number of flux gather invocations") 271 self.lift_counter = EventCounter("n_lift", 272 "Number of flux lift invocations") 273 self.mass_counter = EventCounter("n_mass_op", 274 "Number of mass operator applications") 275 self.diff_counter = EventCounter("n_diff", 276 "Number of differentiation operator applications") 277 278 self.gather_flop_counter = EventCounter("n_flops_gather", 279 "Number of floating point operations in gather") 280 self.lift_flop_counter = EventCounter("n_flops_lift", 281 "Number of floating point operations in lift") 282 self.mass_flop_counter = EventCounter("n_flops_mass", 283 "Number of floating point operations in mass operator") 284 self.diff_flop_counter = EventCounter("n_flops_diff", 285 "Number of floating point operations in diff operator") 286 self.vector_math_flop_counter = EventCounter("n_flops_vector_math", 287 "Number of floating point operations in vector math") 288 289 self.interpolant_counter = EventCounter("n_interp", 290 "Number of interpolant evaluations") 291 292 self.interpolant_timer = IntervalTimer("t_interp", 293 "Time spent evaluating interpolants") 294 295 for op in self.create_op_timers(): 296 mgr.add_quantity(op) 297 298 mgr.add_quantity(self.gather_counter) 299 mgr.add_quantity(self.lift_counter) 300 mgr.add_quantity(self.mass_counter) 301 mgr.add_quantity(self.diff_counter) 302 303 mgr.add_quantity(self.gather_flop_counter) 304 mgr.add_quantity(self.lift_flop_counter) 305 mgr.add_quantity(self.mass_flop_counter) 306 mgr.add_quantity(self.diff_flop_counter) 307 mgr.add_quantity(self.vector_math_flop_counter) 308 309 mgr.add_quantity(self.interpolant_counter) 310 mgr.add_quantity(self.interpolant_timer) 311 312 from pytools.log import time_and_count_function 313 self.interpolate_volume_function = \ 314 time_and_count_function( 315 self.interpolate_volume_function, 316 self.interpolant_timer, 317 self.interpolant_counter) 318 319 self.interpolate_boundary_function = \ 320 time_and_count_function( 321 self.interpolate_boundary_function, 322 self.interpolant_timer, 323 self.interpolant_counter) 324 325 from pytools import single_valued 326 try: 327 order = single_valued(eg.local_discretization.order 328 for eg in self.element_groups) 329 except ValueError: 330 pass 331 else: 332 mgr.set_constant("dg_order", order) 333 334 mgr.set_constant("default_type", self.default_scalar_type.__name__) 335 mgr.set_constant("element_count", len(self.mesh.elements)) 336 mgr.set_constant("node_count", len(self.nodes)) 337 338 for f in self.all_debug_flags(): 339 mgr.set_constant("debug_%s" % f, f in self.debug) 340 341 self.instrumented = True

342 343 # initialization ----------------------------------------------------------

344 - def _build_element_groups_and_nodes(self, local_discretization):

345 from hedge._internal import UniformElementRanges 346 347 eg = _ElementGroup() 348 eg.members = self.mesh.elements 349 eg.member_nrs = numpy.fromiter((el.id for el in eg.members), dtype=numpy.uint32) 350 eg.local_discretization = ldis = local_discretization 351 eg.ranges = UniformElementRanges( 352 0, 353 len(ldis.unit_nodes()), 354 len(self.mesh.elements)) 355 356 nodes_per_el = ldis.node_count() 357 # mem layout: 358 # [....element....][...element...] 359 # | \ 360 # [node.] 361 # | | | 362 # x y z 363 364 # nodes should not be a multi-d array: this would break once 365 # p-adaptivity is implemented 366 self.nodes = numpy.empty( 367 (len(self.mesh.elements)*nodes_per_el, self.dimensions), 368 dtype=float, order="C") 369 370 unit_nodes = numpy.empty( (nodes_per_el, self.dimensions), 371 dtype=float, order="C") 372 373 for i_node, node in enumerate(ldis.unit_nodes()): 374 unit_nodes[i_node] = node 375 376 from hedge._internal import map_element_nodes 377 378 for el in self.mesh.elements: 379 map_element_nodes( 380 self.nodes, 381 el.id*nodes_per_el*self.dimensions, 382 el.map, 383 unit_nodes, 384 self.dimensions) 385 386 self.group_map = [(eg, i) for i in range(len(self.mesh.elements))] 387 self.element_groups = [eg]

388

389 - def _calculate_local_matrices(self):

390 for eg in self.element_groups: 391 ldis = eg.local_discretization 392 393 mmat = eg.mass_matrix = ldis.mass_matrix() 394 immat = eg.inverse_mass_matrix = ldis.inverse_mass_matrix() 395 dmats = eg.differentiation_matrices = \ 396 ldis.differentiation_matrices() 397 smats = eg.stiffness_matrices = \ 398 [numpy.dot(mmat, d) for d in dmats] 399 smats = eg.stiffness_t_matrices = \ 400 [numpy.dot(d.T, mmat.T) for d in dmats] 401 eg.minv_st = \ 402 [numpy.dot(numpy.dot(immat,d.T), mmat) for d in dmats] 403 404 eg.jacobians = numpy.array([ 405 abs(el.map.jacobian()) 406 for el in eg.members]) 407 eg.inverse_jacobians = numpy.array([ 408 abs(el.inverse_map.jacobian()) 409 for el in eg.members]) 410 411 eg.diff_coefficients = numpy.array([ 412 [ 413 [ 414 el.inverse_map 415 .matrix[loc_coord, glob_coord] 416 for el in eg.members 417 ] 418 for loc_coord in range(ldis.dimensions) 419 ] 420 for glob_coord in range(ldis.dimensions) 421 ]) 422 423 eg.stiffness_coefficients = numpy.array([ 424 [ 425 [ 426 abs(el.map.jacobian())*el.inverse_map 427 .matrix[loc_coord, glob_coord] 428 for el in eg.members 429 ] 430 for loc_coord in range(ldis.dimensions) 431 ] 432 for glob_coord in range(ldis.dimensions) 433 ])

434

435 - def _set_flux_face_data(self, f, ldis, (el, fi)):

436 f.face_jacobian = el.face_jacobians[fi] 437 f.element_id = el.id 438 f.face_id = fi 439 f.order = ldis.order 440 f.normal = el.face_normals[fi] 441 442 # This crude approximation is shamelessly stolen from sledge. 443 # There's an important caveat, however (which took me the better 444 # part of a week to figure out): 445 # h on both sides of an interface must be the same, otherwise 446 # the penalty term will behave very oddly. 447 # This unification happens below. 448 f.h = abs(el.map.jacobian()/f.face_jacobian)

449

450 - def _build_interior_face_groups(self):

451 from hedge._internal import FacePair 452 from hedge.element import FaceVertexMismatch 453 454 fg = _FaceGroup(double_sided=True, 455 debug="ilist_generation" in self.debug) 456 457 all_ldis_l = [] 458 all_ldis_n = [] 459 460 debug_node_perm = "node_permutation" in self.debug 461 462 # find and match node indices along faces 463 for i, (local_face, neigh_face) in enumerate(self.mesh.interfaces): 464 e_l, fi_l = local_face 465 e_n, fi_n = neigh_face 466 467 eslice_l, ldis_l = self.find_el_data(e_l.id) 468 eslice_n, ldis_n = self.find_el_data(e_n.id) 469 470 all_ldis_l.append(ldis_l) 471 all_ldis_n.append(ldis_n) 472 473 vertices_l = e_l.faces[fi_l] 474 vertices_n = e_n.faces[fi_n] 475 476 findices_l = ldis_l.face_indices()[fi_l] 477 findices_n = ldis_n.face_indices()[fi_n] 478 479 try: 480 findices_shuffle_op_n = \ 481 ldis_l.get_face_index_shuffle_to_match( 482 vertices_l, vertices_n) 483 484 if debug_node_perm and ldis_l.has_facial_nodes and ldis_n.has_facial_nodes: 485 findices_shuffled_n = findices_shuffle_op_n(findices_n) 486 487 for i, j in zip(findices_l, findices_shuffled_n): 488 dist = self.nodes[eslice_l.start+i]-self.nodes[eslice_n.start+j] 489 assert la.norm(dist) < 1e-14 490 491 except FaceVertexMismatch: 492 # this happens if vertices_l is not a permutation of vertices_n. 493 # periodicity is the only reason why that would be so. 494 495 vertices_n, axis = self.mesh.periodic_opposite_faces[vertices_n] 496 497 findices_shuffle_op_n = \ 498 ldis_l.get_face_index_shuffle_to_match(vertices_l, vertices_n) 499 500 if debug_node_perm and ldis_l.has_facial_nodes and ldis_n.has_facial_nodes: 501 findices_shuffled_n = findices_shuffle_op_n(findices_n) 502 503 for i, j in zip(findices_l, findices_shuffled_n): 504 dist = self.nodes[eslice_l.start+i]-self.nodes[eslice_n.start+j] 505 dist[axis] = 0 506 assert la.norm(dist) < 1e-14 507 508 # create and fill the face pair 509 fp = FacePair() 510 511 fp.loc.el_base_index = eslice_l.start 512 fp.opp.el_base_index = eslice_n.start 513 514 fp.loc.face_index_list_number = fg.register_face_index_list( 515 identifier=fi_l, 516 generator=lambda: findices_l) 517 fp.opp.face_index_list_number = fg.register_face_index_list( 518 identifier=(fi_n, findices_shuffle_op_n), 519 generator=lambda : findices_shuffle_op_n(findices_n)) 520 from pytools import get_write_to_map_from_permutation 521 fp.opp_native_write_map = fg.register_face_index_list( 522 identifier=(fi_n, findices_shuffle_op_n, "wtm"), 523 generator=lambda : 524 get_write_to_map_from_permutation( 525 findices_shuffle_op_n(findices_n), findices_n)) 526 527 self._set_flux_face_data(fp.loc, ldis_l, local_face) 528 self._set_flux_face_data(fp.opp, ldis_n, neigh_face) 529 530 # unify h across the faces 531 fp.loc.h = fp.opp.h = max(fp.loc.h, fp.opp.h) 532 533 assert len(fp.__dict__) == 0 534 assert len(fp.loc.__dict__) == 0 535 assert len(fp.opp.__dict__) == 0 536 537 fg.face_pairs.append(fp) 538 539 if len(fg.face_pairs): 540 from pytools import single_valued 541 ldis_l = single_valued(all_ldis_l) 542 ldis_n = single_valued(all_ldis_n) 543 544 fg.commit(self, ldis_l, ldis_n) 545 546 self.face_groups = [fg] 547 else: 548 self.face_groups = []

549

550 - def boundary_nonempty(self, tag):

551 return bool(self.mesh.tag_to_boundary.get(tag, []))

552 553 @memoize_method

554 - def get_boundary(self, tag):

555 """Get a _Boundary instance for a given `tag'. 556 557 If there is no boundary tagged with `tag', an empty _Boundary instance 558 is returned. Asking for a nonexistant boundary is not an error. 559 (Otherwise get_boundary would unnecessarily become non-local when run 560 in parallel.) 561 """ 562 from hedge._internal import FacePair 563 564 nodes = [] 565 face_ranges = {} 566 vol_indices = [] 567 face_group = _FaceGroup(double_sided=False, 568 debug="ilist_generation" in self.debug) 569 ldis = None # if this boundary is empty, we might as well have no ldis 570 el_face_to_face_group_and_face_pair = {} 571 572 for ef in self.mesh.tag_to_boundary.get(tag, []): 573 el, face_nr = ef 574 575 el_slice, ldis = self.find_el_data(el.id) 576 face_indices = ldis.face_indices()[face_nr] 577 578 f_start = len(nodes) 579 nodes += [self.nodes[el_slice.start+i] for i in face_indices] 580 face_ranges[ef] = (f_start, len(nodes)) 581 vol_indices.extend(el_slice.start+i for i in face_indices) 582 583 # create the face pair 584 fp = FacePair() 585 fp.loc.el_base_index = el_slice.start 586 fp.opp.el_base_index = f_start 587 fp.loc.face_index_list_number = face_group.register_face_index_list( 588 identifier=face_nr, 589 generator=lambda: face_indices) 590 fp.opp.face_index_list_number = face_group.register_face_index_list( 591 identifier=(), 592 generator=lambda: tuple(xrange(len(face_indices)))) 593 self._set_flux_face_data(fp.loc, ldis, ef) 594 assert len(fp.__dict__) == 0 595 assert len(fp.loc.__dict__) == 0 596 assert len(fp.opp.__dict__) == 0 597 598 face_group.face_pairs.append(fp) 599 600 # and make it possible to find it later 601 el_face_to_face_group_and_face_pair[ef] = \ 602 face_group, len(face_group.face_pairs)-1 603 604 face_group.commit(self, ldis, ldis) 605 606 bdry = _Boundary( 607 nodes=numpy.array(nodes), 608 ranges=face_ranges, 609 vol_indices=numpy.asarray(vol_indices, dtype=numpy.intp), 610 face_groups=[face_group], 611 el_face_to_face_group_and_face_pair= 612 el_face_to_face_group_and_face_pair) 613 614 return bdry

615 616 # vector construction -----------------------------------------------------

617 - def __len__(self):

618 """Return the number of nodes in this discretization.""" 619 return len(self.nodes)

620

621 - def len_boundary(self, tag):

622 return len(self.get_boundary(tag).nodes)

623

624 - def get_kind(self, field):

625 return "numpy"

626 627 compute_kind = "numpy" 628

629 - def convert_dtype(self, field, dtype):

630 from hedge.tools import with_object_array_or_scalar 631 if dtype is not None: 632 return with_object_array_or_scalar(lambda f: f.astype(dtype), field) 633 else: 634 return field

635

636 - def convert_volume(self, field, kind):

637 orig_kind = self.get_kind(field) 638 639 if orig_kind != "numpy": 640 raise ValueError, "unable to perform kind conversion: %s -> %s" % ( 641 orig_kind, kind) 642 643 return field

644

645 - def convert_boundary(self, field, tag, kind):

646 orig_kind = self.get_kind(field) 647 648 if orig_kind != "numpy": 649 raise ValueError, "unable to perform kind conversion: %s -> %s" % ( 650 orig_kind, kind) 651 652 return field

653

654 - def convert_boundary_async(self, field, tag, kind, read_map=None):

655 from hedge.tools import ImmediateFuture 656 657 if read_map is not None: 658 from hedge.tools import log_shape 659 ls = log_shape(field) 660 if field.dtype == object or ls == (): 661 from hedge.tools import with_object_array_or_scalar 662 field = with_object_array_or_scalar( 663 lambda f: f[read_map], field) 664 else: 665 field = numpy.asarray( 666 numpy.take(field, read_map, axis=len(ls)), 667 order="C") 668 669 return ImmediateFuture( 670 self.convert_boundary(field, tag, kind))

671

672 - def volume_empty(self, shape=(), dtype=None, kind="numpy"):

673 if kind != "numpy": 674 raise ValueError, "invalid vector kind requested" 675 676 if dtype is None: 677 dtype = self.default_scalar_type 678 return numpy.empty(shape+(len(self.nodes),), dtype)

679

680 - def volume_zeros(self, shape=(), dtype=None, kind="numpy"):

681 if kind != "numpy": 682 raise ValueError, "invalid vector kind requested" 683 684 if dtype is None: 685 dtype = self.default_scalar_type 686 return numpy.zeros(shape+(len(self.nodes),), dtype)

687

688 - def interpolate_volume_function(self, f, dtype=None, kind=None):

689 if kind is None: 690 kind = self.compute_kind 691 692 try: 693 # are we interpolating many fields at once? 694 shape = f.shape 695 except AttributeError: 696 # no, just one 697 shape = () 698 699 slice_pfx = (slice(None),)*len(shape) 700 out = self.volume_empty(shape, dtype, kind="numpy") 701 for eg in self.element_groups: 702 for el, el_slice in zip(eg.members, eg.ranges): 703 for point_nr in xrange(el_slice.start, el_slice.stop): 704 out[slice_pfx + (point_nr,)] = \ 705 f(self.nodes[point_nr], el) 706 return self.convert_volume(out, kind=kind)

707

708 - def boundary_empty(self, tag, shape=(), dtype=None, kind="numpy"):

709 if kind not in ["numpy", "numpy-mpi-recv"]: 710 raise ValueError, "invalid vector kind requested" 711 712 if dtype is None: 713 dtype = self.default_scalar_type 714 return numpy.empty(shape+(len(self.get_boundary(tag).nodes),), dtype)

715

716 - def boundary_zeros(self, tag, shape=(), dtype=None, kind="numpy"):

717 if kind not in ["numpy", "numpy-mpi-recv"]: 718 raise ValueError, "invalid vector kind requested" 719 if dtype is None: 720 dtype = self.default_scalar_type 721 722 return numpy.zeros(shape+(len(self.get_boundary(tag).nodes),), dtype)

723

724 - def interpolate_boundary_function(self, f, tag, dtype=None, kind=None):

725 if kind is None: 726 kind = self.compute_kind 727 728 try: 729 # are we interpolating many fields at once? 730 shape = f.shape 731 except AttributeError: 732 # no, just one 733 shape = () 734 735 out = self.boundary_zeros(tag, shape, dtype, kind="numpy") 736 slice_pfx = (slice(None),)*len(shape) 737 for point_nr, x in enumerate(self.get_boundary(tag).nodes): 738 out[slice_pfx + (point_nr,)] = f(x, None) # FIXME 739 740 return self.convert_boundary(out, tag, kind)

741 742 @memoize_method

743 - def boundary_normals(self, tag, dtype=None, kind=None):

744 if kind is None: 745 kind = self.compute_kind 746 747 result = self.boundary_zeros(shape=(self.dimensions,), tag=tag, dtype=dtype, 748 kind="numpy") 749 for fg in self.get_boundary(tag).face_groups: 750 for face_pair in fg.face_pairs: 751 oeb = face_pair.opp.el_base_index 752 opp_index_list = fg.index_lists[face_pair.opp.face_index_list_number] 753 for i in opp_index_list: 754 result[:,oeb+i] = face_pair.loc.normal 755 756 return self.convert_boundary(result, tag, kind)

757

758 - def volumize_boundary_field(self, bfield, tag, kind=None):

759 if kind is None: 760 kind = self.compute_kind 761 762 if kind != "numpy": 763 raise ValueError("invalid target vector kind in volumize_boundary_field") 764 765 bdry = self.get_boundary(tag) 766 767 def f(subfld): 768 result = self.volume_zeros(dtype=bfield.dtype, kind="numpy") 769 result[bdry.vol_indices] = subfld 770 return result

771 772 from hedge.tools import with_object_array_or_scalar 773 return with_object_array_or_scalar(f, bfield)

774

775 - def boundarize_volume_field(self, field, tag, kind=None):

776 if kind is None: 777 kind = self.compute_kind 778 779 if kind != "numpy": 780 raise ValueError("invalid target vector kind in boundarize_volume_field") 781 782 bdry = self.get_boundary(tag) 783 784 from hedge.tools import log_shape, is_obj_array 785 ls = log_shape(field) 786 787 if is_obj_array(field): 788 if len(field) == 0: 789 return numpy.zeros(()) 790 791 result = self.boundary_empty(tag, shape=ls, dtype=field[0].dtype) 792 from pytools import indices_in_shape 793 for i in indices_in_shape(ls): 794 result[i] = field[i][bdry.vol_indices] 795 796 return result 797 else: 798 return field[tuple(slice(None) for i in range(len(ls))) + (bdry.vol_indices,)]

799

800 - def boundarize_volume_field_async(self, field, tag, kind=None):

801 from hedge.tools import ImmediateFuture 802 return ImmediateFuture( 803 self.boundarize_volume_field(field, tag, kind))

804

805 - def prepare_from_neighbor_map(self, indices):

806 return numpy.array(indices, dtype=numpy.intp)

807 808 # scalar reduction --------------------------------------------------------

809 - def nodewise_dot_product(self, a, b):

810 return numpy.dot(a, b)

811 812 @memoize_method

813 - def _mass_ones(self):

814 from hedge.optemplate import MassOperator 815 return MassOperator().apply(self, ones_on_volume(self))

816 817 @memoize_method

818 - def mesh_volume(self):

819 return self.integral(ones_on_volume(self))

820

821 - def integral(self, volume_vector):

822 from hedge.tools import log_shape 823 824 ls = log_shape(volume_vector) 825 if ls == (): 826 if isinstance(volume_vector, (int, float, complex)): 827 # accept scalars as volume_vector 828 empty = self.volume_empty(dtype=type(volume_vector)) 829 empty.fill(volume_vector) 830 volume_vector = empty 831 832 return self.nodewise_dot_product( 833 self._mass_ones(), volume_vector) 834 else: 835 result = numpy.zeros(shape=ls, dtype=float) 836 837 from pytools import indices_in_shape 838 for i in indices_in_shape(ls): 839 vvi = volume_vector[i] 840 if isinstance(vvi, (int, float)) and vvi == 0: 841 result[i] = 0 842 else: 843 result[i] = self.nodewise_dot_product( 844 self._mass_ones(), volume_vector[i]) 845 846 return result

847 848 @memoize_method

849 - def _compiled_mass_operator(self):

850 from hedge.optemplate import MassOperator, Field 851 mass_op_func = self.compile(MassOperator() * Field("f")) 852 return lambda f: mass_op_func(f=f)

853

854 - def norm(self, volume_vector, p=2):

855 if p == numpy.Inf: 856 return numpy.abs(volume_vector).max() 857 else: 858 from hedge.tools import log_shape 859 860 if p != 2: 861 volume_vector = numpy.abs(volume_vector)**(p/2) 862 863 return self.inner_product( 864 volume_vector, 865 volume_vector)**(1/p)

866

867 - def inner_product(self, a, b):

868 mass_op = self._compiled_mass_operator() 869 870 from hedge.tools import log_shape 871 ls = log_shape(a) 872 assert log_shape(b) == ls 873 if ls == (): 874 return float(self.nodewise_dot_product( 875 a, mass_op(b))) 876 else: 877 assert len(ls) == 1 878 return float(sum( 879 self.nodewise_dot_product( 880 sub_a, mass_op(sub_b)) 881 for sub_a, sub_b in zip(a,b)))

882 883 # element data retrieval --------------------------------------------------

884 - def find_el_range(self, el_id):

885 group, idx = self.group_map[el_id] 886 return group.ranges[idx]

887

888 - def find_el_discretization(self, el_id):

889 return self.group_map[el_id][0].local_discretization

890

891 - def find_el_data(self, el_id):

892 group, idx = self.group_map[el_id] 893 return group.ranges[idx], group.local_discretization

894

895 - def find_element(self, idx):

896 for i, (start, stop) in enumerate(self.element_group): 897 if start <= idx < stop: 898 return i 899 raise ValueError, "not a valid dof index"

900 901 # misc stuff --------------------------------------------------------------

902 - def dt_non_geometric_factor(self):

903 distinct_ldis = set(eg.local_discretization for eg in self.element_groups) 904 return min(ldis.dt_non_geometric_factor() 905 for ldis in distinct_ldis)

906

907 - def dt_geometric_factor(self):

908 return min(min(eg.local_discretization.dt_geometric_factor( 909 [self.mesh.points[i] for i in el.vertex_indices], el) 910 for el in eg.members) 911 for eg in self.element_groups)

912

913 - def dt_factor(self, max_system_ev, stepper_class=None, *stepper_args):

914 u"""Calculate the largest stable timestep, given a time stepper 915 `stepper_class`. If none is given, RK4 is assumed. 916 """ 917 918 # Calculating the correct timestep Δt for a DG scheme using the RK4 919 # method is described in: "Nodal DG Methods, Algorithm, Analysis and 920 # Applications" by J.S. Hesthaven & T. Warburton, p. 93, "Discrete 921 # stability and timestep choise". The implementation of timestep 922 # calculation here is based upon this chapter. 923 # 924 # For a spatially continuous problem, the timestep can be calculated by 925 # the following relation: 926 # 927 # max|λop| * Δt = C_TimeStepper, 928 # 929 # where max|λop| is the maximum eigenvalue of the operator and 930 # C_TimeStepper represents the maximum size of the stability region of 931 # the timestepper along the imaginary axis. 932 # 933 # For a DG-discretized problem another factor has to be added: 934 # 935 # fDG = fNG * fG, 936 # 937 # fNG: non geometric factor fG: geometric factor 938 # 939 # The discrete relation is: max|λop| * Δt = fDG * C_Timestepper 940 # 941 # Since the LocalDiscretization.dt_non_geometric_factor() and 942 # LocalDiscretization.dt_geometric_factor() implicitly scale their 943 # results for an RK4 time stepper, fDG includes already C_RK4 such as 944 # fDG becomes fDG_RK4 and the relation is: 945 # 946 # max|λop| * Δt = fDG_RK4 947 # 948 # As this is only sufficient for the use of RK4 timestepper but not for 949 # any other implemented approache (e.g. Adams-Bashforth) additional 950 # information about the size of the stability region is required to be 951 # added into the relation. 952 # 953 # Unifying the relation with the size of the RK4 stability region and 954 # multiplying it with the size of the specific timestepper stability 955 # region brings out the correct relation: 956 # 957 # max|λop| * Δt = fDG_RK4 / C_RK4 * C_TimeStepper 958 # 959 # C_TimeStepper gets calculated by a bisection method for every kind of 960 # timestepper. 961 962 963 rk4_dt = 1/max_system_ev \ 964 * self.dt_non_geometric_factor() \ 965 * self.dt_geometric_factor() 966 967 from hedge.timestep import RK4TimeStepper 968 if stepper_class is None or stepper_class == RK4TimeStepper: 969 return rk4_dt 970 else: 971 assert isinstance(stepper_class, type) 972 973 from hedge.timestep.stability import \ 974 calculate_fudged_stability_region 975 976 return rk4_dt \ 977 * calculate_fudged_stability_region( 978 stepper_class, *stepper_args) \ 979 / calculate_fudged_stability_region(RK4TimeStepper)

980

981 - def get_point_evaluator(self, point):

982 for eg in self.element_groups: 983 for el, rng in zip(eg.members, eg.ranges): 984 if el.contains_point(point): 985 ldis = eg.local_discretization 986 basis_values = numpy.array([ 987 phi(el.inverse_map(point)) 988 for phi in ldis.basis_functions()]) 989 vdm_t = ldis.vandermonde().T 990 return _PointEvaluator( 991 discr=self, 992 el_range=rng, 993 interp_coeff=la.solve(vdm_t, basis_values)) 994 995 raise RuntimeError, "point %s not found" % point

996 997 # op template execution ---------------------------------------------------

998 - def compile(self, optemplate, post_bind_mapper=lambda x: x):

999 ex = self.executor_class(self, optemplate, post_bind_mapper) 1000 1001 if "dump_dataflow_graph" in self.debug: 1002 ex.code.dump_dataflow_graph() 1003 1004 if self.instrumented: 1005 ex.instrument() 1006 return ex

1007

1008 - def add_function(self, name, func):

1009 self.exec_functions[name] = func

1010

1011 1012 1013 1014 # random utilities ------------------------------------------------------------ 1015 -class SymmetryMap(object):

1016 """A symmetry map on global DG functions. 1017 1018 Suppose that the L{Mesh} on which a L{Discretization} is defined has 1019 is mapped onto itself by a nontrivial symmetry map M{f(.)}. Then 1020 this class allows you to carry out this map on vectors representing 1021 functions on this L{Discretization}. 1022 """

1023 - def __init__(self, discr, sym_map, element_map, threshold=1e-13):

1024 self.discretization = discr 1025 1026 complete_el_map = {} 1027 for i, j in element_map.iteritems(): 1028 complete_el_map[i] = j 1029 complete_el_map[j] = i 1030 1031 self.map = {} 1032 1033 for eg in discr.element_groups: 1034 for el, el_slice in zip(eg.members, eg.ranges): 1035 mapped_i_el = complete_el_map[el.id] 1036 mapped_slice = discr.find_el_range(mapped_i_el) 1037 for i_pt in range(el_slice.start, el_slice.stop): 1038 pt = discr.nodes[i_pt] 1039 mapped_pt = sym_map(pt) 1040 for m_i_pt in range(mapped_slice.start, mapped_slice.stop): 1041 if la.norm(discr.nodes[m_i_pt] - mapped_pt) < threshold: 1042 self.map[i_pt] = m_i_pt 1043 break 1044 1045 if i_pt not in self.map: 1046 for m_i_pt in range(mapped_slice.start, mapped_slice.stop): 1047 print la.norm_2(discr.nodes[m_i_pt] - mapped_pt) 1048 raise RuntimeError, "no symmetry match found"

1049

1050 - def __call__(self, vec):

1051 result = self.discretization.volume_zeros() 1052 for i, mapped_i in self.map.iteritems(): 1053 result[mapped_i] = vec[i] 1054 return result

1055

1056 1057 1058 1059 -def generate_random_constant_on_elements(discr):

1060 result = discr.volume_zeros() 1061 import random 1062 for eg in discr.element_groups: 1063 for e_start, e_end in eg.ranges: 1064 result[e_start:e_end] = random.random() 1065 return result

1066

1067 1068 1069 1070 -def ones_on_boundary(discr, tag):

1071 result = discr.volume_zeros() 1072 1073 try: 1074 faces = discr.mesh.tag_to_boundary[tag] 1075 except KeyError: 1076 pass 1077 else: 1078 for face in faces: 1079 el, fl = face 1080 1081 (el_start, el_end), ldis = discr.find_el_data(el.id) 1082 fl_indices = ldis.face_indices()[fl] 1083 1084 for i in fl_indices: 1085 result[el_start+i] = 1 1086 1087 return result

1088

1089 1090 1091 1092 -def ones_on_volume(discr):

1093 result = discr.volume_empty() 1094 result.fill(1) 1095 return result

1096

1097 1098 1099 1100 # projection between different discretizations -------------------------------- 1101 -class Projector:

1102 - def __init__(self, from_discr, to_discr):

1103 self.from_discr = from_discr 1104 self.to_discr = to_discr 1105 1106 self.interp_matrices = [] 1107 for from_eg, to_eg in zip( 1108 from_discr.element_groups, to_discr.element_groups): 1109 from_ldis = from_eg.local_discretization 1110 to_ldis = to_eg.local_discretization 1111 1112 from_count = from_ldis.node_count() 1113 to_count = to_ldis.node_count() 1114 1115 # check that the two element groups have the same members 1116 for from_el, to_el in zip(from_eg.members, to_eg.members): 1117 assert from_el is to_el 1118 1119 from hedge.tools import permutation_matrix 1120 1121 # assemble the from->to mode permutation matrix, guided by 1122 # mode identifiers 1123 if to_count > from_count: 1124 to_node_ids_to_idx = dict( 1125 (nid, i) for i, nid in 1126 enumerate(to_ldis.generate_mode_identifiers())) 1127 1128 to_indices = [ 1129 to_node_ids_to_idx[from_nid] 1130 for from_nid in from_ldis.generate_mode_identifiers() 1131 ] 1132 1133 pmat = permutation_matrix( 1134 to_indices=to_indices, 1135 h=to_count, w=from_count) 1136 else: 1137 from_node_ids_to_idx = dict( 1138 (nid, i) for i, nid in 1139 enumerate(from_ldis.generate_mode_identifiers())) 1140 1141 from_indices = [ 1142 from_node_ids_to_idx[to_nid] 1143 for to_nid in to_ldis.generate_mode_identifiers() 1144 ] 1145 1146 pmat = permutation_matrix( 1147 from_indices=from_indices, 1148 h=to_count, w=from_count) 1149 1150 # build interpolation matrix 1151 from_matrix = from_ldis.vandermonde() 1152 to_matrix = to_ldis.vandermonde() 1153 1154 from hedge.tools import leftsolve 1155 from numpy import dot 1156 self.interp_matrices.append( 1157 numpy.asarray( 1158 leftsolve(from_matrix, dot(to_matrix, pmat)), 1159 order="C"))

1160

1161 - def __call__(self, from_vec):

1162 from hedge._internal import perform_elwise_operator 1163 from hedge.tools import log_shape 1164 1165 ls = log_shape(from_vec) 1166 result = self.to_discr.volume_zeros(ls, kind="numpy") 1167 1168 from pytools import indices_in_shape 1169 for i in indices_in_shape(ls): 1170 for from_eg, to_eg, imat in zip( 1171 self.from_discr.element_groups, 1172 self.to_discr.element_groups, 1173 self.interp_matrices): 1174 perform_elwise_operator( 1175 from_eg.ranges, to_eg.ranges, 1176 imat, from_vec[i], result[i]) 1177 1178 return result

1179

1180 1181 1182 1183 # filter ---------------------------------------------------------------------- 1184 -class ExponentialFilterResponseFunction:

1185 """A typical exponential-falloff mode response filter function. 1186 1187 See description in Section 5.6.1 of Hesthaven/Warburton. 1188 """

1189 - def __init__(self, min_amplification=0.1, order=6):

1190 """Construct the filter function. 1191 1192 The amplification factor of the lowest-order (constant) mode is always 1. 1193 1194 @arg min_amplification: The amplification factor applied to the highest mode. 1195 @arg order: The order of the filter. This controls how fast (or slowly) the 1196 C{min_amplification} is reached. 1197 """ 1198 from math import log 1199 self.alpha = -log(min_amplification) 1200 self.order = order

1201

1202 - def __call__(self, mode_idx, ldis):

1203 eta = sum(mode_idx)/ldis.order 1204 1205 from math import exp 1206 return exp(-self.alpha * eta**self.order)

1207

1208 1209 1210 1211 -class Filter:

1212 - def __init__(self, discr, mode_response_func):

1213 """Construct a filter. 1214 1215 @arg discr: The L{Discretization} for which the filter is to be 1216 constructed. 1217 @arg mode_response_func: A function mapping 1218 C{(mode_tuple, local_discretization)} to a float indicating the 1219 factor by which this mode is to be multiplied after filtering. 1220 """ 1221 self.discr = discr 1222 1223 self.filter_map = {} 1224 1225 for eg in discr.element_groups: 1226 ldis = eg.local_discretization 1227 1228 node_count = ldis.node_count() 1229 1230 filter_coeffs = [mode_response_func(mid, ldis) 1231 for mid in ldis.generate_mode_identifiers()] 1232 1233 # build filter matrix 1234 vdm = ldis.vandermonde() 1235 from hedge.tools import leftsolve 1236 from numpy import dot 1237 mat = numpy.asarray( 1238 leftsolve(vdm, 1239 dot(vdm, numpy.diag(filter_coeffs))), 1240 order="C") 1241 self.filter_map[eg] = mat

1242

1243 - def __call__(self, vec):

1244 from hedge.tools import log_shape 1245 1246 ls = log_shape(vec) 1247 result = self.discr.volume_zeros(ls) 1248 1249 from pytools import indices_in_shape 1250 for i in indices_in_shape(ls): 1251 from hedge._internal import perform_elwise_operator 1252 for eg in self.discr.element_groups: 1253 perform_elwise_operator(eg.ranges, eg.ranges, 1254 self.filter_map[eg], vec[i], result[i]) 1255 1256 return result

1257

1258 - def get_filter_matrix(self, el_group):

1259 return self.filter_map[el_group]

1260

Source Code for Module hedge.discretization