diff --git a/examples/wave/wave-min.py b/examples/wave/wave-min.py index 42f3041eebe8f7b3603028c8edabfa7a4cbc9af1..6a5a6bf907abef0c7a12c55cefe5de8881af8ce8 100644 --- a/examples/wave/wave-min.py +++ b/examples/wave/wave-min.py @@ -28,6 +28,7 @@ THE SOFTWARE. import numpy as np import pyopencl as cl from grudge.shortcuts import make_discretization, set_up_rk4 +from grudge import sym, bind def main(write_output=True): @@ -43,7 +44,6 @@ def main(write_output=True): source_width = 0.05 source_omega = 3 - from grudge import sym sym_x = sym.nodes(2) sym_source_center_dist = sym_x - source_center sym_sin = sym.CFunction("sin") @@ -71,6 +71,14 @@ def main(write_output=True): # FIXME #dt = op.estimate_rk4_timestep(discr, fields=fields) + op.check_bc_coverage(mesh) + + print(sym.pretty(op.sym_operator())) + bound_op = bind(discr, op.sym_operator()) + + def rhs(t, w): + return bound_op(t=t, w=w) + dt = 0.001 dt_stepper = set_up_rk4(dt, fields, rhs) diff --git a/grudge/__init__.py b/grudge/__init__.py index 38a84cd2918b14d6805bcf1506b0d1d8d40bc145..1b765e560e7d9f5f86f47689c52aacaf37d97fbd 100644 --- a/grudge/__init__.py +++ b/grudge/__init__.py @@ -22,5 +22,7 @@ OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN THE SOFTWARE. """ -import grudge.symbolic.primitives as sym # noqa +import grudge.sym # noqa import grudge.symbolic.flux.primitives as sym_flux # noqa + +from grudge.execution import bind # noqa diff --git a/grudge/execution.py b/grudge/execution.py new file mode 100644 index 0000000000000000000000000000000000000000..0452eb3200be5d3b4856de9cb71306d2c4b29db8 --- /dev/null +++ b/grudge/execution.py @@ -0,0 +1,674 @@ +from __future__ import division, absolute_import + +__copyright__ = "Copyright (C) 2015 Andreas Kloeckner" + +__license__ = """ +Permission is hereby granted, free of charge, to any person obtaining a copy +of this software and associated documentation files (the "Software"), to deal +in the Software without restriction, including without limitation the rights +to use, copy, modify, merge, publish, distribute, sublicense, and/or sell +copies of the Software, and to permit persons to whom the Software is +furnished to do so, subject to the following conditions: + +The above copyright notice and this permission notice shall be included in +all copies or substantial portions of the Software. + +THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR +IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY, +FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL THE +AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER +LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM, +OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN +THE SOFTWARE. +""" + +import numpy as np +import grudge.symbolic.mappers as mappers + +import logging +logger = logging.getLogger(__name__) + + +# {{{ exec mapper + +class ExecutionMapper(mappers.Evaluator, + mappers.BoundOpMapperMixin, + mappers.LocalOpReducerMixin): + def __init__(self, context, executor): + super(ExecutionMapper, self).__init__(context) + self.discr = executor.discr + self.executor = executor + + # {{{ expression mappings ------------------------------------------------- + + def map_ones(self, expr): + # FIXME + if expr.quadrature_tag is not None: + raise NotImplementedError("ones on quad. grids") + + result = self.discr.volume_empty(kind=self.discr.compute_kind) + result.fill(1) + return result + + def map_node_coordinate_component(self, expr): + # FIXME + if expr.quadrature_tag is not None: + raise NotImplementedError("node coordinate components on quad. grids") + # FIXME: Data transfer and strided CPU index every time, ugh. + return self.discr.convert_volume( + # Yes, that .copy() is necessary because much of grudge + # doesn't check for striding. Ugh^2. + self.discr.nodes[:, expr.axis].copy(), + kind=self.discr.compute_kind) + + def map_normal_component(self, expr): + if expr.quadrature_tag is not None: + raise NotImplementedError("normal components on quad. grids") + return self.discr.boundary_normals(expr.boundary_tag)[expr.axis] + + def map_boundarize(self, op, field_expr): + return self.discr.boundarize_volume_field( + self.rec(field_expr), tag=op.tag, + kind=self.discr.compute_kind) + + def map_scalar_parameter(self, expr): + return self.context[expr.name] + + def map_jacobian(self, expr): + return self.discr.volume_jacobians(expr.quadrature_tag) + + def map_forward_metric_derivative(self, expr): + return (self.discr.forward_metric_derivatives(expr.quadrature_tag) + [expr.xyz_axis][expr.rst_axis]) + + def map_inverse_metric_derivative(self, expr): + return (self.discr.inverse_metric_derivatives(expr.quadrature_tag) + [expr.xyz_axis][expr.rst_axis]) + + def map_call(self, expr): + from pymbolic.primitives import Variable + assert isinstance(expr.function, Variable) + func_name = expr.function.name + + try: + func = self.discr.exec_functions[func_name] + except KeyError: + func = getattr(np, expr.function.name) + + return func(*[self.rec(p) for p in expr.parameters]) + + def map_nodal_sum(self, op, field_expr): + return np.sum(self.rec(field_expr)) + + def map_nodal_max(self, op, field_expr): + return np.max(self.rec(field_expr)) + + def map_nodal_min(self, op, field_expr): + return np.min(self.rec(field_expr)) + + def map_if_positive(self, expr): + crit = self.rec(expr.criterion) + bool_crit = crit > 0 + then = self.rec(expr.then) + else_ = self.rec(expr.else_) + + true_indices = np.nonzero(bool_crit) + false_indices = np.nonzero(~bool_crit) + + result = np.empty_like(crit) + + if isinstance(then, np.ndarray): + then = then[true_indices] + if isinstance(else_, np.ndarray): + else_ = else_[false_indices] + + result[true_indices] = then + result[false_indices] = else_ + return result + + def map_if(self, expr): + bool_crit = self.rec(expr.condition) + then = self.rec(expr.then) + else_ = self.rec(expr.else_) + + true_indices = np.nonzero(bool_crit) + false_indices = np.nonzero(~bool_crit) + + result = self.discr.volume_empty( + kind=self.discr.compute_kind) + + if isinstance(then, np.ndarray): + then = then[true_indices] + if isinstance(else_, np.ndarray): + else_ = else_[false_indices] + + result[true_indices] = then + result[false_indices] = else_ + return result + + def map_ref_diff_base(self, op, field_expr): + raise NotImplementedError( + "differentiation should be happening in batched form") + + def map_elementwise_linear(self, op, field_expr): + field = self.rec(field_expr) + + from grudge.tools import is_zero + if is_zero(field): + return 0 + + out = self.discr.volume_zeros() + self.executor.do_elementwise_linear(op, field, out) + return out + + def map_ref_quad_mass(self, op, field_expr): + field = self.rec(field_expr) + + from grudge.tools import is_zero + if is_zero(field): + return 0 + + qtag = op.quadrature_tag + + from grudge._internal import perform_elwise_operator + + out = self.discr.volume_zeros() + for eg in self.discr.element_groups: + eg_quad_info = eg.quadrature_info[qtag] + + perform_elwise_operator(eg_quad_info.ranges, eg.ranges, + eg_quad_info.ldis_quad_info.mass_matrix(), + field, out) + + return out + + def map_quad_grid_upsampler(self, op, field_expr): + field = self.rec(field_expr) + + from grudge.tools import is_zero + if is_zero(field): + return 0 + + qtag = op.quadrature_tag + + from grudge._internal import perform_elwise_operator + quad_info = self.discr.get_quadrature_info(qtag) + + out = np.zeros(quad_info.node_count, field.dtype) + for eg in self.discr.element_groups: + eg_quad_info = eg.quadrature_info[qtag] + + perform_elwise_operator(eg.ranges, eg_quad_info.ranges, + eg_quad_info.ldis_quad_info.volume_up_interpolation_matrix(), + field, out) + + return out + + def map_quad_int_faces_grid_upsampler(self, op, field_expr): + field = self.rec(field_expr) + + from grudge.tools import is_zero + if is_zero(field): + return 0 + + qtag = op.quadrature_tag + + from grudge._internal import perform_elwise_operator + quad_info = self.discr.get_quadrature_info(qtag) + + out = np.zeros(quad_info.int_faces_node_count, field.dtype) + for eg in self.discr.element_groups: + eg_quad_info = eg.quadrature_info[qtag] + + perform_elwise_operator(eg.ranges, eg_quad_info.el_faces_ranges, + eg_quad_info.ldis_quad_info.volume_to_face_up_interpolation_matrix() + .copy(), + field, out) + + return out + + def map_quad_bdry_grid_upsampler(self, op, field_expr): + field = self.rec(field_expr) + + from grudge.tools import is_zero + if is_zero(field): + return 0 + + bdry = self.discr.get_boundary(op.boundary_tag) + bdry_q_info = bdry.get_quadrature_info(op.quadrature_tag) + + out = np.zeros(bdry_q_info.node_count, field.dtype) + + from grudge._internal import perform_elwise_operator + for fg, from_ranges, to_ranges, ldis_quad_info in zip( + bdry.face_groups, + bdry.fg_ranges, + bdry_q_info.fg_ranges, + bdry_q_info.fg_ldis_quad_infos): + perform_elwise_operator(from_ranges, to_ranges, + ldis_quad_info.face_up_interpolation_matrix(), + field, out) + + return out + + def map_elementwise_max(self, op, field_expr): + from grudge._internal import perform_elwise_max + field = self.rec(field_expr) + + out = self.discr.volume_zeros(dtype=field.dtype) + for eg in self.discr.element_groups: + perform_elwise_max(eg.ranges, field, out) + + return out + + # }}} + + # {{{ code execution functions -------------------------------------------- + def exec_assign(self, insn): + return [(name, self.rec(expr)) + for name, expr in zip(insn.names, insn.exprs)], [] + + def exec_vector_expr_assign(self, insn): + if self.discr.instrumented: + def stats_callback(n, vec_expr): + self.discr.vector_math_flop_counter.add(n*insn.flop_count()) + return self.discr.vector_math_timer + else: + stats_callback = None + + if insn.flop_count() == 0: + return [(name, self(expr)) + for name, expr in zip(insn.names, insn.exprs)], [] + else: + compiled = insn.compiled(self.executor) + return zip(compiled.result_names(), + compiled(self, stats_callback)), [] + + def exec_flux_batch_assign(self, insn): + from pymbolic.primitives import is_zero + + class ZeroSpec: + pass + + class BoundaryZeros(ZeroSpec): + pass + + class VolumeZeros(ZeroSpec): + pass + + def eval_arg(arg_spec): + arg_expr, is_int = arg_spec + arg = self.rec(arg_expr) + if is_zero(arg): + if insn.is_boundary and not is_int: + return BoundaryZeros() + else: + return VolumeZeros() + else: + return arg + + args = [eval_arg(arg_expr) + for arg_expr in insn.flux_var_info.arg_specs] + + from pytools import common_dtype + max_dtype = common_dtype( + [a.dtype for a in args if not isinstance(a, ZeroSpec)], + self.discr.default_scalar_type) + + def cast_arg(arg): + if isinstance(arg, BoundaryZeros): + return self.discr.boundary_zeros( + insn.repr_op.boundary_tag, dtype=max_dtype) + elif isinstance(arg, VolumeZeros): + return self.discr.volume_zeros( + dtype=max_dtype) + elif isinstance(arg, np.ndarray): + return np.asarray(arg, dtype=max_dtype) + else: + return arg + + args = [cast_arg(arg) for arg in args] + + if insn.quadrature_tag is None: + if insn.is_boundary: + face_groups = self.discr.get_boundary(insn.repr_op.boundary_tag)\ + .face_groups + else: + face_groups = self.discr.face_groups + else: + if insn.is_boundary: + face_groups = self.discr.get_boundary(insn.repr_op.boundary_tag)\ + .get_quadrature_info(insn.quadrature_tag).face_groups + else: + face_groups = self.discr.get_quadrature_info(insn.quadrature_tag) \ + .face_groups + + result = [] + + for fg in face_groups: + # grab module + module = insn.get_module(self.discr, max_dtype) + func = module.gather_flux + + # set up argument structure + arg_struct = module.ArgStruct() + for arg_name, arg in zip(insn.flux_var_info.arg_names, args): + setattr(arg_struct, arg_name, arg) + for arg_num, scalar_arg_expr in enumerate( + insn.flux_var_info.scalar_parameters): + setattr(arg_struct, + "_scalar_arg_%d" % arg_num, + self.rec(scalar_arg_expr)) + + fof_shape = (fg.face_count*fg.face_length()*fg.element_count(),) + all_fluxes_on_faces = [ + np.zeros(fof_shape, dtype=max_dtype) + for f in insn.expressions] + for i, fof in enumerate(all_fluxes_on_faces): + setattr(arg_struct, "flux%d_on_faces" % i, fof) + + # make sure everything ended up in Boost.Python attributes + # (i.e. empty __dict__) + assert not arg_struct.__dict__, arg_struct.__dict__.keys() + + # perform gather + func(fg, arg_struct) + + # do lift, produce output + for name, flux_bdg, fluxes_on_faces in zip(insn.names, insn.expressions, + all_fluxes_on_faces): + + if insn.quadrature_tag is None: + if flux_bdg.op.is_lift: + mat = fg.ldis_loc.lifting_matrix() + scaling = fg.local_el_inverse_jacobians + else: + mat = fg.ldis_loc.multi_face_mass_matrix() + scaling = None + else: + assert not flux_bdg.op.is_lift + mat = fg.ldis_loc_quad_info.multi_face_mass_matrix() + scaling = None + + out = self.discr.volume_zeros(dtype=fluxes_on_faces.dtype) + self.executor.lift_flux(fg, mat, scaling, fluxes_on_faces, out) + + if self.discr.instrumented: + from grudge.tools import lift_flops + + # correct for quadrature, too. + self.discr.lift_flop_counter.add(lift_flops(fg)) + + result.append((name, out)) + + if not face_groups: + # No face groups? Still assign context variables. + for name, flux_bdg in zip(insn.names, insn.expressions): + result.append((name, self.discr.volume_zeros())) + + return result, [] + + def exec_diff_batch_assign(self, insn): + rst_diff = self.executor.diff(insn.operators, self.rec(insn.field)) + + return [(name, diff) for name, diff in zip(insn.names, rst_diff)], [] + + exec_quad_diff_batch_assign = exec_diff_batch_assign + + # }}} + +# }}} + + +# {{{ executor + +class Executor(object): + def __init__(self, discr, code, debug_flags): + self.discr = discr + self.code = code + self.elwise_linear_cache = {} + self.debug_flags = debug_flags + + if "dump_op_code" in debug_flags: + from grudge.tools import open_unique_debug_file + open_unique_debug_file("op-code", ".txt").write( + str(self.code)) + + def bench_diff(f): + test_field = discr.volume_zeros() + from grudge.optemplate import ReferenceDifferentiationOperator + from time import time + + start = time() + f([ReferenceDifferentiationOperator(i) + for i in range(discr.dimensions)], test_field) + return time() - start + + def bench_lift(f): + if len(discr.face_groups) == 0: + return 0 + + fg = discr.face_groups[0] + out = discr.volume_zeros() + from time import time + + fof_shape = (fg.face_count*fg.face_length()*fg.element_count(),) + fof = np.zeros(fof_shape, dtype=self.discr.default_scalar_type) + + start = time() + f(fg, fg.ldis_loc.lifting_matrix(), + fg.local_el_inverse_jacobians, fof, out) + return time() - start + + def pick_faster_func(benchmark, choices, attempts=3): + from pytools import argmin2 + return argmin2( + (f, min(benchmark(f) for i in range(attempts))) + for f in choices) + + from grudge.backends.jit.diff import JitDifferentiator + self.diff = pick_faster_func(bench_diff, + [self.diff_builtin, JitDifferentiator(discr)]) + from grudge.backends.jit.lift import JitLifter + self.lift_flux = pick_faster_func(bench_lift, + [self.lift_flux, JitLifter(discr)]) + + def instrument(self): + discr = self.discr + assert discr.instrumented + + from pytools.log import time_and_count_function + from grudge.tools import time_count_flop + + from grudge.tools import diff_rst_flops, mass_flops + + if discr.quad_min_degrees: + from warnings import warn + warn("flop counts for quadrature may be wrong") + + self.diff_rst = \ + time_count_flop( + self.diff_rst, + discr.diff_timer, + discr.diff_counter, + discr.diff_flop_counter, + diff_rst_flops(discr)) + + self.do_elementwise_linear = \ + time_count_flop( + self.do_elementwise_linear, + discr.el_local_timer, + discr.el_local_counter, + discr.el_local_flop_counter, + mass_flops(discr)) + + self.lift_flux = \ + time_and_count_function( + self.lift_flux, + discr.lift_timer, + discr.lift_counter) + + def lift_flux(self, fgroup, matrix, scaling, field, out): + from grudge._internal import lift_flux + from pytools import to_uncomplex_dtype + lift_flux(fgroup, + matrix.astype(to_uncomplex_dtype(field.dtype)), + scaling, field, out) + + def diff_rst(self, op, field): + result = self.discr.volume_zeros(dtype=field.dtype) + + from grudge._internal import perform_elwise_operator + for eg in self.discr.element_groups: + perform_elwise_operator(op.preimage_ranges(eg), eg.ranges, + op.matrices(eg)[op.rst_axis].astype(field.dtype), + field, result) + + return result + + def diff_builtin(self, operators, field): + """For the batch of reference differentiation operators in + *operators*, return the local corresponding derivatives of + *field*. + """ + + return [self.diff_rst(op, field) for op in operators] + + def do_elementwise_linear(self, op, field, out): + for eg in self.discr.element_groups: + try: + matrix, coeffs = self.elwise_linear_cache[eg, op, field.dtype] + except KeyError: + matrix = np.asarray(op.matrix(eg), dtype=field.dtype) + coeffs = op.coefficients(eg) + self.elwise_linear_cache[eg, op, field.dtype] = matrix, coeffs + + from grudge._internal import ( + perform_elwise_scaled_operator, + perform_elwise_operator) + + if coeffs is None: + perform_elwise_operator(eg.ranges, eg.ranges, + matrix, field, out) + else: + perform_elwise_scaled_operator(eg.ranges, eg.ranges, + coeffs, matrix, field, out) + + def __call__(self, **context): + return self.code.execute( + self.discr.exec_mapper_class(context, self)) + +# }}} + + +# {{{ process_sym_operator function + +def process_sym_operator(sym_operator, post_bind_mapper=None, + dumper=lambda name, sym_operator: None, mesh=None, + type_hints={}): + + from grudge.symbolic.mappers import ( + OperatorBinder, CommutativeConstantFoldingMapper, + EmptyFluxKiller, InverseMassContractor, DerivativeJoiner, + ErrorChecker, OperatorSpecializer, GlobalToReferenceMapper) + from grudge.symbolic.mappers.bc_to_flux import BCToFluxRewriter + from grudge.symbolic.mappers.type_inference import TypeInferrer + + dumper("before-bind", sym_operator) + sym_operator = OperatorBinder()(sym_operator) + + ErrorChecker(mesh)(sym_operator) + + if post_bind_mapper is not None: + dumper("before-postbind", sym_operator) + sym_operator = post_bind_mapper(sym_operator) + + if mesh is not None: + dumper("before-empty-flux-killer", sym_operator) + sym_operator = EmptyFluxKiller(mesh)(sym_operator) + + dumper("before-cfold", sym_operator) + sym_operator = CommutativeConstantFoldingMapper()(sym_operator) + + dumper("before-bc2flux", sym_operator) + sym_operator = BCToFluxRewriter()(sym_operator) + + # Ordering restriction: + # + # - Must run constant fold before first type inference pass, because zeros, + # while allowed, violate typing constraints (because they can't be assigned + # a unique type), and need to be killed before the type inferrer sees them. + # + # - Must run BC-to-flux before first type inferrer run so that zeros in + # flux arguments can be removed. + + dumper("before-specializer", sym_operator) + sym_operator = OperatorSpecializer( + TypeInferrer()(sym_operator, type_hints) + )(sym_operator) + + # Ordering restriction: + # + # - Must run OperatorSpecializer before performing the GlobalToReferenceMapper, + # because otherwise it won't differentiate the type of grids (node or quadrature + # grids) that the operators will apply on. + + assert mesh is not None + dumper("before-global-to-reference", sym_operator) + sym_operator = GlobalToReferenceMapper(mesh.ambient_dim)(sym_operator) + + # Ordering restriction: + # + # - Must specialize quadrature operators before performing inverse mass + # contraction, because there are no inverse-mass-contracted variants of the + # quadrature operators. + + dumper("before-imass", sym_operator) + sym_operator = InverseMassContractor()(sym_operator) + + dumper("before-cfold-2", sym_operator) + sym_operator = CommutativeConstantFoldingMapper()(sym_operator) + + dumper("before-derivative-join", sym_operator) + sym_operator = DerivativeJoiner()(sym_operator) + + dumper("process-finished", sym_operator) + + return sym_operator + +# }}} + + +def bind(discr, sym_operator, post_bind_mapper=lambda x: x, type_hints={}, + debug_flags=set()): + # from grudge.symbolic.mappers import QuadratureUpsamplerRemover + # sym_operator = QuadratureUpsamplerRemover(self.quad_min_degrees)( + # sym_operator) + + stage = [0] + + def dump_optemplate(name, sym_operator): + if "dump_optemplate_stages" in debug_flags: + from grudge.tools import open_unique_debug_file + from grudge.optemplate import pretty + open_unique_debug_file("%02d-%s" % (stage[0], name), ".txt").write( + pretty(sym_operator)) + stage[0] += 1 + + sym_operator = process_sym_operator(sym_operator, + post_bind_mapper=post_bind_mapper, + dumper=dump_optemplate, + mesh=discr.mesh, + type_hints=type_hints) + + from grudge.symbolic.compiler import OperatorCompiler + code = OperatorCompiler(discr)(sym_operator, type_hints) + + ex = Executor(discr, code, type_hints) + + if "dump_dataflow_graph" in debug_flags: + ex.code.dump_dataflow_graph() + + return ex + +# vim: foldmethod=marker diff --git a/grudge/models/wave.py b/grudge/models/wave.py index c7671259e3cc02a9b7c9938e2b39df96618928d8..f016699347aaa1fe3655e0b34170a2b7ce254592 100644 --- a/grudge/models/wave.py +++ b/grudge/models/wave.py @@ -85,7 +85,7 @@ class StrongWaveOperator(HyperbolicOperator): w = sym_flux.FluxVectorPlaceholder(1+dim) u = w[0] v = w[1:] - normal = sym_flux.make_normal(dim) + normal = sym_flux.normal(dim) flux_weak = join_fields( np.dot(v.avg, normal), @@ -147,7 +147,7 @@ class StrongWaveOperator(HyperbolicOperator): ) # entire operator ----------------------------------------------------- - nabla = sym.make_nabla(d) + nabla = sym.nabla(d) flux_op = sym.get_flux_operator(self.flux()) result = ( @@ -167,20 +167,13 @@ class StrongWaveOperator(HyperbolicOperator): return result - def bind(self, discr): - from grudge.mesh import check_bc_coverage - check_bc_coverage(discr.mesh, [ + def check_bc_coverage(self, mesh): + from meshmode.mesh import check_bc_coverage + check_bc_coverage(mesh, [ self.dirichlet_tag, self.neumann_tag, self.radiation_tag]) - compiled_sym_operator = discr.compile(self.op_template()) - - def rhs(t, w, **extra_context): - return compiled_sym_operator(t=t, w=w, **extra_context) - - return rhs - def max_eigenvalue(self, t, fields=None, discr=None): return abs(self.c) @@ -232,14 +225,12 @@ class VariableVelocityStrongWaveOperator(HyperbolicOperator): # {{{ flux ---------------------------------------------------------------- def flux(self): - from grudge.flux import FluxVectorPlaceholder, make_normal - dim = self.dimensions - w = FluxVectorPlaceholder(2+dim) + w = sym_flux.FluxVectorPlaceholder(2+dim) c = w[0] u = w[1] v = w[2:] - normal = make_normal(dim) + normal = sym_flux.normal(dim) flux = self.time_sign*1/2*join_fields( c.ext * np.dot(v.ext, normal) @@ -272,18 +263,9 @@ class VariableVelocityStrongWaveOperator(HyperbolicOperator): return do def sym_operator(self, with_sensor=False): - from grudge.symbolic import \ - Field, \ - make_sym_vector, \ - BoundaryPair, \ - get_flux_operator, \ - make_nabla, \ - InverseMassOperator, \ - BoundarizeOperator - d = self.dimensions - w = make_sym_vector("w", d+1) + w = sym.make_sym_vector("w", d+1) u = w[0] v = w[1:] @@ -292,26 +274,25 @@ class VariableVelocityStrongWaveOperator(HyperbolicOperator): # {{{ boundary conditions # Dirichlet - dir_c = BoundarizeOperator(self.dirichlet_tag) * self.c - dir_u = BoundarizeOperator(self.dirichlet_tag) * u - dir_v = BoundarizeOperator(self.dirichlet_tag) * v + dir_c = sym.BoundarizeOperator(self.dirichlet_tag) * self.c + dir_u = sym.BoundarizeOperator(self.dirichlet_tag) * u + dir_v = sym.BoundarizeOperator(self.dirichlet_tag) * v dir_bc = join_fields(dir_c, -dir_u, dir_v) # Neumann - neu_c = BoundarizeOperator(self.neumann_tag) * self.c - neu_u = BoundarizeOperator(self.neumann_tag) * u - neu_v = BoundarizeOperator(self.neumann_tag) * v + neu_c = sym.BoundarizeOperator(self.neumann_tag) * self.c + neu_u = sym.BoundarizeOperator(self.neumann_tag) * u + neu_v = sym.BoundarizeOperator(self.neumann_tag) * v neu_bc = join_fields(neu_c, neu_u, -neu_v) # Radiation - from grudge.symbolic import make_normal - rad_normal = make_normal(self.radiation_tag, d) + rad_normal = sym.make_normal(self.radiation_tag, d) - rad_c = BoundarizeOperator(self.radiation_tag) * self.c - rad_u = BoundarizeOperator(self.radiation_tag) * u - rad_v = BoundarizeOperator(self.radiation_tag) * v + rad_c = sym.BoundarizeOperator(self.radiation_tag) * self.c + rad_u = sym.BoundarizeOperator(self.radiation_tag) * u + rad_v = sym.BoundarizeOperator(self.radiation_tag) * v rad_bc = join_fields( rad_c, @@ -333,7 +314,7 @@ class VariableVelocityStrongWaveOperator(HyperbolicOperator): diffusion_coeff = self.diffusion_coeff if with_sensor: - diffusion_coeff += Field("sensor") + diffusion_coeff += sym.Field("sensor") from grudge.second_order import SecondDerivativeTarget @@ -360,8 +341,8 @@ class VariableVelocityStrongWaveOperator(HyperbolicOperator): # }}} # entire operator ----------------------------------------------------- - nabla = make_nabla(d) - flux_op = get_flux_operator(self.flux()) + nabla = sym.nabla(d) + flux_op = sym.get_flux_operator(self.flux()) return ( - join_fields( @@ -372,32 +353,20 @@ class VariableVelocityStrongWaveOperator(HyperbolicOperator): make_diffusion, v) ) + - InverseMassOperator() * ( + sym.InverseMassOperator() * ( flux_op(flux_w) - + flux_op(BoundaryPair(flux_w, dir_bc, self.dirichlet_tag)) - + flux_op(BoundaryPair(flux_w, neu_bc, self.neumann_tag)) - + flux_op(BoundaryPair(flux_w, rad_bc, self.radiation_tag)) + + flux_op(sym.BoundaryPair(flux_w, dir_bc, self.dirichlet_tag)) + + flux_op(sym.BoundaryPair(flux_w, neu_bc, self.neumann_tag)) + + flux_op(sym.BoundaryPair(flux_w, rad_bc, self.radiation_tag)) )) - def bind(self, discr, sensor=None): - from grudge.mesh import check_bc_coverage - check_bc_coverage(discr.mesh, [ + def check_bc_coverage(self, mesh): + from meshmode.mesh import check_bc_coverage + check_bc_coverage(mesh, [ self.dirichlet_tag, self.neumann_tag, self.radiation_tag]) - compiled_sym_operator = discr.compile(self.op_template( - with_sensor=sensor is not None)) - - def rhs(t, w): - kwargs = {} - if sensor is not None: - kwargs["sensor"] = sensor(t, w) - - return compiled_sym_operator(t=t, w=w, **kwargs) - - return rhs - def max_eigenvalue_expr(self): import grudge.symbolic as sym return sym.NodalMax()(sym.CFunction("fabs")(self.c)) diff --git a/grudge/sym.py b/grudge/sym.py new file mode 100644 index 0000000000000000000000000000000000000000..bf9c7de615c82c11a5a6d91fa9ed872341d0d32c --- /dev/null +++ b/grudge/sym.py @@ -0,0 +1,28 @@ +from __future__ import division, absolute_import + +__copyright__ = "Copyright (C) 2015 Andreas Kloeckner" + +__license__ = """ +Permission is hereby granted, free of charge, to any person obtaining a copy +of this software and associated documentation files (the "Software"), to deal +in the Software without restriction, including without limitation the rights +to use, copy, modify, merge, publish, distribute, sublicense, and/or sell +copies of the Software, and to permit persons to whom the Software is +furnished to do so, subject to the following conditions: + +The above copyright notice and this permission notice shall be included in +all copies or substantial portions of the Software. + +THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR +IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY, +FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL THE +AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER +LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM, +OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN +THE SOFTWARE. +""" + +from grudge.symbolic.primitives import * # noqa +from grudge.symbolic.operators import * # noqa + +from grudge.symbolic.tools import pretty # noqa diff --git a/grudge/symbolic/compiler.py b/grudge/symbolic/compiler.py index c1375e220a1b3254e82f779784af8d9a1c8857ab..e8b0abbf690bdc075c862524b0072f846ed03545 100644 --- a/grudge/symbolic/compiler.py +++ b/grudge/symbolic/compiler.py @@ -1,13 +1,8 @@ """Compiler to turn operator expression tree into (imperative) bytecode.""" -from __future__ import division -from __future__ import absolute_import -from __future__ import print_function -import six -from six.moves import zip -from functools import reduce +from __future__ import division, absolute_import, print_function -__copyright__ = "Copyright (C) 2008 Andreas Kloeckner" +__copyright__ = "Copyright (C) 2008-15 Andreas Kloeckner" __license__ = """ Permission is hereby granted, free of charge, to any person obtaining a copy @@ -30,8 +25,11 @@ THE SOFTWARE. """ +import six +from six.moves import zip, reduce from pytools import Record, memoize_method -from grudge.symbolic import IdentityMapper +from grudge import sym +from grudge.symbolic.mappers import IdentityMapper # {{{ instructions @@ -549,7 +547,7 @@ class Code(object): # {{{ compiler -class OperatorCompilerBase(IdentityMapper): +class OperatorCompiler(IdentityMapper): class FluxRecord(Record): __slots__ = ["flux_expr", "dependencies", "repr_op"] @@ -601,8 +599,7 @@ class OperatorCompilerBase(IdentityMapper): # {{{ top-level driver ---------------------------------------------------- def __call__(self, expr, type_hints={}): # Put the result expressions into variables as well. - from grudge.symbolic import make_common_subexpression as cse - expr = cse(expr, "_result") + expr = sym.make_common_subexpression(expr, "_result") from grudge.symbolic.mappers.type_inference import TypeInferrer self.typedict = TypeInferrer()(expr, type_hints) @@ -833,10 +830,10 @@ class OperatorCompilerBase(IdentityMapper): try: return self.expr_to_var[expr] except KeyError: - from grudge.tools import is_field_equal + from pytools.obj_array import obj_array_equal all_flux_xchgs = [fe for fe in self.flux_exchange_ops - if is_field_equal(fe.arg_fields, expr.arg_fields)] + if obj_array_equal(fe.arg_fields, expr.arg_fields)] assert len(all_flux_xchgs) > 0 diff --git a/grudge/symbolic/flux/primitives.py b/grudge/symbolic/flux/primitives.py index 11eb9218492868dbc21495e82c21da256e8356e4..e2629fbe43f12bed850792a46f37d0b14b619f50 100644 --- a/grudge/symbolic/flux/primitives.py +++ b/grudge/symbolic/flux/primitives.py @@ -165,7 +165,7 @@ def norm(v): return numpy.dot(v, v)**0.5 -def make_normal(dimensions): +def normal(dimensions): return numpy.array([Normal(i) for i in range(dimensions)], dtype=object) diff --git a/grudge/symbolic/mappers/__init__.py b/grudge/symbolic/mappers/__init__.py index 5c88cdbbc92647213c0ab045be82e9ce062f7249..3806f9a7d7309b110bc926e5c3614fe4c89686bc 100644 --- a/grudge/symbolic/mappers/__init__.py +++ b/grudge/symbolic/mappers/__init__.py @@ -1,4 +1,4 @@ -"""Operator template mappers.""" +"""Mappers to transform symbolic operators.""" from __future__ import division @@ -35,6 +35,10 @@ import pymbolic.mapper.constant_folder import pymbolic.mapper.flop_counter from pymbolic.mapper import CSECachingMapperMixin +from grudge import sym +from grudge import sym_flux +import grudge.symbolic.operators as op + # {{{ mixins @@ -316,10 +320,9 @@ class OperatorBinder(CSECachingMapperMixin, IdentityMapper): return expr from pymbolic.primitives import flattened_product, Product - from grudge.optemplate import Operator, OperatorBinding first = expr.children[0] - if isinstance(first, Operator): + if isinstance(first, op.Operator): prod = flattened_product(expr.children[1:]) if isinstance(prod, Product) and len(prod.children) > 1: from warnings import warn @@ -327,7 +330,7 @@ class OperatorBinder(CSECachingMapperMixin, IdentityMapper): "operand in a product is ambiguous - " "use the parenthesized form instead." % first) - return OperatorBinding(first, self.rec(prod)) + return sym.OperatorBinding(first, self.rec(prod)) else: return first * self.rec(flattened_product(expr.children[1:])) @@ -338,7 +341,7 @@ class OperatorBinder(CSECachingMapperMixin, IdentityMapper): class OperatorSpecializer(CSECachingMapperMixin, IdentityMapper): """Guided by a typedict obtained through type inference (i.e. by - :class:`grudge.optemplate.mappers.type_inference.TypeInferrrer`), + :class:`grudge.symbolic.mappers.type_inference.TypeInferrrer`), substitutes more specialized operators for generic ones. For example, if type inference has determined that a differentiation @@ -350,7 +353,7 @@ class OperatorSpecializer(CSECachingMapperMixin, IdentityMapper): def __init__(self, typedict): """ :param typedict: generated by - :class:`grudge.optemplate.mappers.type_inference.TypeInferrer`. + :class:`grudge.symbolic.mappers.type_inference.TypeInferrer`. """ self.typedict = typedict @@ -358,25 +361,9 @@ class OperatorSpecializer(CSECachingMapperMixin, IdentityMapper): IdentityMapper.map_common_subexpression def map_operator_binding(self, expr): - from grudge.optemplate.primitives import BoundaryPair - - from grudge.optemplate.operators import ( - MassOperator, - QuadratureMassOperator, - ReferenceMassOperator, - ReferenceQuadratureMassOperator, - - StiffnessTOperator, - QuadratureStiffnessTOperator, - ReferenceStiffnessTOperator, - ReferenceQuadratureStiffnessTOperator, + from grudge.symbolic.primitives import BoundaryPair - QuadratureGridUpsampler, QuadratureBoundaryGridUpsampler, - FluxOperatorBase, FluxOperator, QuadratureFluxOperator, - BoundaryFluxOperator, QuadratureBoundaryFluxOperator, - BoundarizeOperator) - - from grudge.optemplate.mappers.type_inference import ( + from grudge.symbolic.mappers.type_inference import ( type_info, QuadratureRepresentation) # {{{ figure out field type @@ -399,30 +386,31 @@ class OperatorSpecializer(CSECachingMapperMixin, IdentityMapper): QuadratureRepresentation) # }}} - # Based on where this is run in the optemplate processing pipeline, - # we may encounter both reference and non-reference operators. + # Based on where this is run in the symbolic operator processing + # pipeline, we may encounter both reference and non-reference + # operators. # {{{ elementwise operators - if isinstance(expr.op, MassOperator) and has_quad_operand: - return QuadratureMassOperator( + if isinstance(expr.op, op.MassOperator) and has_quad_operand: + return op.QuadratureMassOperator( field_repr_tag.quadrature_tag)(self.rec(expr.field)) - elif isinstance(expr.op, ReferenceMassOperator) and has_quad_operand: - return ReferenceQuadratureMassOperator( + elif isinstance(expr.op, op.ReferenceMassOperator) and has_quad_operand: + return op.ReferenceQuadratureMassOperator( field_repr_tag.quadrature_tag)(self.rec(expr.field)) - elif (isinstance(expr.op, StiffnessTOperator) and has_quad_operand): - return QuadratureStiffnessTOperator( + elif (isinstance(expr.op, op.StiffnessTOperator) and has_quad_operand): + return op.QuadratureStiffnessTOperator( expr.op.xyz_axis, field_repr_tag.quadrature_tag)( self.rec(expr.field)) - elif (isinstance(expr.op, ReferenceStiffnessTOperator) + elif (isinstance(expr.op, op.ReferenceStiffnessTOperator) and has_quad_operand): - return ReferenceQuadratureStiffnessTOperator( + return op.ReferenceQuadratureStiffnessTOperator( expr.op.xyz_axis, field_repr_tag.quadrature_tag)( self.rec(expr.field)) - elif (isinstance(expr.op, QuadratureGridUpsampler) + elif (isinstance(expr.op, op.QuadratureGridUpsampler) and isinstance(field_type, type_info.BoundaryVectorBase)): # potential shortcut: #if (isinstance(expr.field, OperatorBinding) @@ -431,16 +419,16 @@ class OperatorSpecializer(CSECachingMapperMixin, IdentityMapper): #expr.field.op.tag, expr.op.quadrature_tag)( #self.rec(expr.field.field)) - return QuadratureBoundaryGridUpsampler( + return op.QuadratureBoundaryGridUpsampler( expr.op.quadrature_tag, field_type.boundary_tag)(expr.field) # }}} - elif isinstance(expr.op, BoundarizeOperator) and has_quad_operand: + elif isinstance(expr.op, op.BoundarizeOperator) and has_quad_operand: raise TypeError("BoundarizeOperator cannot be applied to " "quadrature-based operands--use QuadUpsample(Boundarize(...))") # {{{ flux operator specialization - elif isinstance(expr.op, FluxOperatorBase): + elif isinstance(expr.op, op.FluxOperatorBase): from pytools.obj_array import with_object_array_or_scalar repr_tag_cell = [None] @@ -473,29 +461,29 @@ class OperatorSpecializer(CSECachingMapperMixin, IdentityMapper): if is_boundary: if is_quad: - return QuadratureBoundaryFluxOperator( + return op.QuadratureBoundaryFluxOperator( flux, quad_tag, bpair.tag)(new_fld) else: - return BoundaryFluxOperator(flux, bpair.tag)(new_fld) + return op.BoundaryFluxOperator(flux, bpair.tag)(new_fld) else: if is_quad: - return QuadratureFluxOperator(flux, quad_tag)(new_fld) + return op.QuadratureFluxOperator(flux, quad_tag)(new_fld) else: - return FluxOperator(flux, expr.op.is_lift)(new_fld) + return op.FluxOperator(flux, expr.op.is_lift)(new_fld) # }}} else: return IdentityMapper.map_operator_binding(self, expr) def map_normal_component(self, expr): - from grudge.optemplate.mappers.type_inference import ( + from grudge.symbolic.mappers.type_inference import ( NodalRepresentation) expr_type = self.typedict[expr] if not isinstance( expr_type.repr_tag, NodalRepresentation): - from grudge.optemplate.primitives import ( + from grudge.symbolic.primitives import ( BoundaryNormalComponent) # for now, parts of this are implemented. @@ -528,27 +516,9 @@ class GlobalToReferenceMapper(CSECachingMapperMixin, IdentityMapper): IdentityMapper.map_common_subexpression def map_operator_binding(self, expr): - from grudge.optemplate.primitives import ( + from grudge.symbolic.primitives import ( Jacobian, InverseMetricDerivative) - from grudge.optemplate.operators import ( - MassOperator, - ReferenceMassOperator, - QuadratureMassOperator, - ReferenceQuadratureMassOperator, - - StiffnessOperator, - - StiffnessTOperator, - ReferenceStiffnessTOperator, - QuadratureStiffnessTOperator, - ReferenceQuadratureStiffnessTOperator, - - InverseMassOperator, ReferenceInverseMassOperator, - DifferentiationOperator, ReferenceDifferentiationOperator, - - MInvSTOperator) - # Global-to-reference is run after operator specialization, so # if we encounter non-quadrature operators here, we know they # must be nodal. @@ -567,43 +537,43 @@ class GlobalToReferenceMapper(CSECachingMapperMixin, IdentityMapper): ref_class(rst_axis, **diff_kwargs)(rec_field) for rst_axis in range(self.dimensions)) - if isinstance(expr.op, MassOperator): - return ReferenceMassOperator()( + if isinstance(expr.op, op.MassOperator): + return op.ReferenceMassOperator()( Jacobian(None) * self.rec(expr.field)) - elif isinstance(expr.op, QuadratureMassOperator): - return ReferenceQuadratureMassOperator( + elif isinstance(expr.op, op.QuadratureMassOperator): + return op.ReferenceQuadratureMassOperator( expr.op.quadrature_tag)( Jacobian(expr.op.quadrature_tag) * self.rec(expr.field)) - elif isinstance(expr.op, InverseMassOperator): - return ReferenceInverseMassOperator()( + elif isinstance(expr.op, op.InverseMassOperator): + return op.ReferenceInverseMassOperator()( 1/Jacobian(None) * self.rec(expr.field)) - elif isinstance(expr.op, StiffnessOperator): - return ReferenceMassOperator()(Jacobian(None) * + elif isinstance(expr.op, op.StiffnessOperator): + return op.ReferenceMassOperator()(Jacobian(None) * self.rec( - DifferentiationOperator(expr.op.xyz_axis)(expr.field))) + op.DifferentiationOperator(expr.op.xyz_axis)(expr.field))) - elif isinstance(expr.op, DifferentiationOperator): + elif isinstance(expr.op, op.DifferentiationOperator): return rewrite_derivative( - ReferenceDifferentiationOperator, + op.ReferenceDifferentiationOperator, expr.field, with_jacobian=False) - elif isinstance(expr.op, StiffnessTOperator): + elif isinstance(expr.op, op.StiffnessTOperator): return rewrite_derivative( - ReferenceStiffnessTOperator, + op.ReferenceStiffnessTOperator, expr.field) - elif isinstance(expr.op, QuadratureStiffnessTOperator): + elif isinstance(expr.op, op.QuadratureStiffnessTOperator): return rewrite_derivative( - ReferenceQuadratureStiffnessTOperator, + op.ReferenceQuadratureStiffnessTOperator, expr.field, quadrature_tag=expr.op.quadrature_tag) - elif isinstance(expr.op, MInvSTOperator): + elif isinstance(expr.op, op.MInvSTOperator): return self.rec( - InverseMassOperator()( - StiffnessTOperator(expr.op.xyz_axis)(expr.field))) + op.InverseMassOperator()( + op.StiffnessTOperator(expr.op.xyz_axis)(expr.field))) else: return IdentityMapper.map_operator_binding(self, expr) @@ -615,8 +585,7 @@ class GlobalToReferenceMapper(CSECachingMapperMixin, IdentityMapper): class StringifyMapper(pymbolic.mapper.stringifier.StringifyMapper): def _format_btag(self, tag): - from grudge.mesh import SYSTEM_TAGS - if tag in SYSTEM_TAGS: + if isinstance(tag, type): return tag.__name__ else: return repr(tag) @@ -626,7 +595,7 @@ class StringifyMapper(pymbolic.mapper.stringifier.StringifyMapper): self, constant_mapper=constant_mapper) if flux_stringify_mapper is None: - from grudge.flux import FluxStringifyMapper + from grudge.symbolic.flux.mappers import FluxStringifyMapper flux_stringify_mapper = FluxStringifyMapper() self.flux_stringify_mapper = flux_stringify_mapper @@ -757,7 +726,7 @@ class StringifyMapper(pymbolic.mapper.stringifier.StringifyMapper): return "ElWMax" def map_boundarize(self, expr, enclosing_prec): - return "Boundarize" % expr.tag + return "Boundarize" % self._format_btag(expr.tag) def map_flux_exchange(self, expr, enclosing_prec): from pymbolic.mapper.stringifier import PREC_NONE @@ -781,10 +750,11 @@ class StringifyMapper(pymbolic.mapper.stringifier.StringifyMapper): def map_normal_component(self, expr, enclosing_prec): if expr.quadrature_tag is None: return ("Normal[%d]" - % (expr.boundary_tag, expr.axis)) + % (self._format_btag(expr.boundary_tag), expr.axis)) else: return ("Q[%s]Normal[%d]" - % (expr.quadrature_tag, expr.boundary_tag, expr.axis)) + % (expr.quadrature_tag, self._format_btag(expr.boundary_tag), + expr.axis)) def map_jacobian(self, expr, enclosing_prec): if expr.quadrature_tag is None: @@ -841,7 +811,7 @@ class PrettyStringifyMapper( self.bc_to_number = {} self.bc_string_list = [] - from grudge.flux import PrettyFluxStringifyMapper + from grudge.symbolic.flux.mappers import PrettyFluxStringifyMapper self.flux_stringify_mapper = PrettyFluxStringifyMapper() def get_flux_number(self, flux): @@ -867,16 +837,15 @@ class PrettyStringifyMapper( self.bc_string_list.append(str_bc) self.bc_to_number[expr] = bc_number - return "BC%d@%s" % (bc_number, expr.tag) + return "BC%d@%s" % (bc_number, self._format_btag(expr.tag)) def map_operator_binding(self, expr, enclosing_prec): - from grudge.optemplate import BoundarizeOperator - if isinstance(expr.op, BoundarizeOperator): + if isinstance(expr.op, op.BoundarizeOperator): from pymbolic.mapper.stringifier import PREC_CALL, PREC_SUM return self.parenthesize_if_needed( "%s@%s" % ( self.rec(expr.field, PREC_CALL), - expr.op.tag), + self._format_btag(expr.op.tag)), enclosing_prec, PREC_SUM) else: return StringifyMapper.map_operator_binding( @@ -935,14 +904,9 @@ class QuadratureUpsamplerRemover(CSECachingMapperMixin, IdentityMapper): IdentityMapper.map_common_subexpression def map_operator_binding(self, expr): - from grudge.optemplate.operators import ( - QuadratureGridUpsampler, - QuadratureInteriorFacesGridUpsampler, - QuadratureBoundaryGridUpsampler) - - if isinstance(expr.op, (QuadratureGridUpsampler, - QuadratureInteriorFacesGridUpsampler, - QuadratureBoundaryGridUpsampler)): + if isinstance(expr.op, (op.QuadratureGridUpsampler, + op.QuadratureInteriorFacesGridUpsampler, + op.QuadratureBoundaryGridUpsampler)): try: min_degree = self.quad_min_degrees[expr.op.quadrature_tag] except KeyError: @@ -984,18 +948,12 @@ class QuadratureDetector(CSECachingMapperMixin, CombineMapper): return self.QuadStatusNotKnown def map_operator_binding(self, expr): - from grudge.optemplate.operators import ( - DiffOperatorBase, FluxOperatorBase, - MassOperatorBase, - QuadratureGridUpsampler, - QuadratureInteriorFacesGridUpsampler) - if isinstance(expr.op, ( - DiffOperatorBase, FluxOperatorBase, - MassOperatorBase)): + op.DiffOperatorBase, op.FluxOperatorBase, + op.MassOperatorBase)): return None - elif isinstance(expr.op, (QuadratureGridUpsampler, - QuadratureInteriorFacesGridUpsampler)): + elif isinstance(expr.op, (op.QuadratureGridUpsampler, + op.QuadratureInteriorFacesGridUpsampler)): return expr.op else: return CombineMapper.map_operator_binding(self, expr) @@ -1016,18 +974,12 @@ class QuadratureUpsamplerChanger(CSECachingMapperMixin, IdentityMapper): IdentityMapper.map_common_subexpression def map_operator_binding(self, expr): - from grudge.optemplate.operators import ( - DiffOperatorBase, FluxOperatorBase, - MassOperatorBase, - QuadratureGridUpsampler, - QuadratureInteriorFacesGridUpsampler) - if isinstance(expr.op, ( - DiffOperatorBase, FluxOperatorBase, - MassOperatorBase)): + op.DiffOperatorBase, op.FluxOperatorBase, + op.MassOperatorBase)): return expr - elif isinstance(expr.op, (QuadratureGridUpsampler, - QuadratureInteriorFacesGridUpsampler)): + elif isinstance(expr.op, (op.QuadratureGridUpsampler, + op.QuadratureInteriorFacesGridUpsampler)): return self.desired_quad_op(expr.field) else: return IdentityMapper.map_operator_binding(self, expr) @@ -1056,11 +1008,8 @@ class CommutativeConstantFoldingMapper( if is_zero(field): return 0 - from grudge.optemplate.operators import FluxOperatorBase - from grudge.optemplate.primitives import BoundaryPair - - if isinstance(expr.op, FluxOperatorBase): - if isinstance(field, BoundaryPair): + if isinstance(expr.op, op.FluxOperatorBase): + if isinstance(field, sym.BoundaryPair): return self.remove_zeros_from_boundary_flux(expr.op, field) else: return self.remove_zeros_from_interior_flux(expr.op, field) @@ -1069,14 +1018,13 @@ class CommutativeConstantFoldingMapper( # {{{ remove zeros from interior flux def remove_zeros_from_interior_flux(self, op, vol_field): - from pytools.obj_array import is_obj_array + from pytools.obj_array import is_obj_array, make_obj_array if not is_obj_array(vol_field): vol_field = [vol_field] - from grudge.flux import FieldComponent subst_map = {} - from grudge.tools import is_zero, make_obj_array + from grudge.tools import is_zero # process volume field new_vol_field = [] @@ -1085,21 +1033,21 @@ class CommutativeConstantFoldingMapper( flux_arg = self.rec(flux_arg) if is_zero(flux_arg): - subst_map[FieldComponent(i, is_interior=True)] = 0 - subst_map[FieldComponent(i, is_interior=False)] = 0 + subst_map[sym_flux.FieldComponent(i, is_interior=True)] = 0 + subst_map[sym_flux.FieldComponent(i, is_interior=False)] = 0 else: new_vol_field.append(flux_arg) - subst_map[FieldComponent(i, is_interior=True)] = \ - FieldComponent(new_idx, is_interior=True) - subst_map[FieldComponent(i, is_interior=False)] = \ - FieldComponent(new_idx, is_interior=False) + subst_map[sym_flux.FieldComponent(i, is_interior=True)] = \ + sym_flux.FieldComponent(new_idx, is_interior=True) + subst_map[sym_flux.FieldComponent(i, is_interior=False)] = \ + sym_flux.FieldComponent(new_idx, is_interior=False) new_idx += 1 # substitute results into flux def sub_flux(expr): return subst_map.get(expr, expr) - from grudge.flux import FluxSubstitutionMapper + from grudge.symbolic.flux.mappers import FluxSubstitutionMapper new_flux = FluxSubstitutionMapper(sub_flux)(op.flux) if is_zero(new_flux): @@ -1114,55 +1062,54 @@ class CommutativeConstantFoldingMapper( def remove_zeros_from_boundary_flux(self, op, bpair): vol_field = bpair.field bdry_field = bpair.bfield - from pytools.obj_array import is_obj_array + from pytools.obj_array import is_obj_array, make_obj_array if not is_obj_array(vol_field): vol_field = [vol_field] if not is_obj_array(bdry_field): bdry_field = [bdry_field] - from grudge.flux import FieldComponent subst_map = {} # process volume field - from grudge.tools import is_zero, make_obj_array + from grudge.tools import is_zero new_vol_field = [] new_idx = 0 for i, flux_arg in enumerate(vol_field): - fc = FieldComponent(i, is_interior=True) + fc = sym_flux.FieldComponent(i, is_interior=True) flux_arg = self.rec(flux_arg) if is_zero(flux_arg): subst_map[fc] = 0 else: new_vol_field.append(flux_arg) - subst_map[fc] = FieldComponent(new_idx, is_interior=True) + subst_map[fc] = sym_flux.FieldComponent(new_idx, is_interior=True) new_idx += 1 # process boundary field new_bdry_field = [] new_idx = 0 for i, flux_arg in enumerate(bdry_field): - fc = FieldComponent(i, is_interior=False) + fc = sym_flux.FieldComponent(i, is_interior=False) flux_arg = self.rec(flux_arg) if is_zero(flux_arg): subst_map[fc] = 0 else: new_bdry_field.append(flux_arg) - subst_map[fc] = FieldComponent(new_idx, is_interior=False) + subst_map[fc] = sym_flux.FieldComponent(new_idx, is_interior=False) new_idx += 1 # substitute results into flux def sub_flux(expr): return subst_map.get(expr, expr) - from grudge.flux import FluxSubstitutionMapper + from grudge.symbolic.flux.mappers import FluxSubstitutionMapper new_flux = FluxSubstitutionMapper(sub_flux)(op.flux) if is_zero(new_flux): return 0 else: - from grudge.optemplate.primitives import BoundaryPair + from grudge.symbolic.primitives import BoundaryPair return type(op)(new_flux, *op.__getinitargs__()[1:])( BoundaryPair( make_obj_array(new_vol_field), @@ -1181,9 +1128,7 @@ class EmptyFluxKiller(CSECachingMapperMixin, IdentityMapper): IdentityMapper.map_common_subexpression def map_operator_binding(self, expr): - from grudge.optemplate import BoundaryFluxOperatorBase - - if (isinstance(expr.op, BoundaryFluxOperatorBase) and + if (isinstance(expr.op, op.BoundaryFluxOperatorBase) and len(self.mesh.tag_to_boundary.get(expr.op.boundary_tag, [])) == 0): return 0 else: @@ -1192,9 +1137,7 @@ class EmptyFluxKiller(CSECachingMapperMixin, IdentityMapper): class _InnerDerivativeJoiner(pymbolic.mapper.RecursiveMapper): def map_operator_binding(self, expr, derivatives): - from grudge.optemplate import DifferentiationOperator - - if isinstance(expr.op, DifferentiationOperator): + if isinstance(expr.op, op.DifferentiationOperator): derivatives.setdefault(expr.op, []).append(expr.field) return 0 else: @@ -1210,7 +1153,7 @@ class _InnerDerivativeJoiner(pymbolic.mapper.RecursiveMapper): self.rec(child, derivatives) for child in expr.children)) def map_product(self, expr, derivatives): - from grudge.optemplate.tools import is_scalar + from grudge.symbolic.tools import is_scalar from pytools import partition scalars, nonscalars = partition(is_scalar, expr.children) @@ -1308,61 +1251,50 @@ class _InnerInverseMassContractor(pymbolic.mapper.RecursiveMapper): def map_constant(self, expr): from grudge.tools import is_zero - from grudge.optemplate import InverseMassOperator, OperatorBinding if is_zero(expr): return 0 else: - return OperatorBinding( - InverseMassOperator(), + return op.OperatorBinding( + op.InverseMassOperator(), self.outer_mass_contractor(expr)) def map_algebraic_leaf(self, expr): - from grudge.optemplate import InverseMassOperator, OperatorBinding - - return OperatorBinding( - InverseMassOperator(), + return op.OperatorBinding( + op.InverseMassOperator(), self.outer_mass_contractor(expr)) def map_operator_binding(self, binding): - from grudge.optemplate import ( - MassOperator, StiffnessOperator, StiffnessTOperator, - DifferentiationOperator, - MInvSTOperator, InverseMassOperator, - FluxOperator, BoundaryFluxOperator) - - if isinstance(binding.op, MassOperator): + if isinstance(binding.op, op.MassOperator): return binding.field - elif isinstance(binding.op, StiffnessOperator): - return DifferentiationOperator(binding.op.xyz_axis)( + elif isinstance(binding.op, op.StiffnessOperator): + return op.DifferentiationOperator(binding.op.xyz_axis)( self.outer_mass_contractor(binding.field)) - elif isinstance(binding.op, StiffnessTOperator): - return MInvSTOperator(binding.op.xyz_axis)( + elif isinstance(binding.op, op.StiffnessTOperator): + return op.MInvSTOperator(binding.op.xyz_axis)( self.outer_mass_contractor(binding.field)) - elif isinstance(binding.op, FluxOperator): + elif isinstance(binding.op, op.FluxOperator): assert not binding.op.is_lift - return FluxOperator(binding.op.flux, is_lift=True)( + return op.FluxOperator(binding.op.flux, is_lift=True)( self.outer_mass_contractor(binding.field)) - elif isinstance(binding.op, BoundaryFluxOperator): + elif isinstance(binding.op, op.BoundaryFluxOperator): assert not binding.op.is_lift - return BoundaryFluxOperator(binding.op.flux, + return op.BoundaryFluxOperator(binding.op.flux, binding.op.boundary_tag, is_lift=True)( self.outer_mass_contractor(binding.field)) else: self.extra_operator_count += 1 - return InverseMassOperator()( + return op.InverseMassOperator()( self.outer_mass_contractor(binding)) def map_sum(self, expr): return expr.__class__(tuple(self.rec(child) for child in expr.children)) def map_product(self, expr): - from grudge.optemplate import InverseMassOperator, ScalarParameter - def is_scalar(expr): - return isinstance(expr, (int, float, complex, ScalarParameter)) + return isinstance(expr, (int, float, complex, sym.ScalarParameter)) from pytools import len_iterable nonscalar_count = len_iterable(ch @@ -1372,7 +1304,7 @@ class _InnerInverseMassContractor(pymbolic.mapper.RecursiveMapper): if nonscalar_count > 1: # too complicated, don't touch it self.extra_operator_count += 1 - return InverseMassOperator()( + return op.InverseMassOperator()( self.outer_mass_contractor(expr)) else: def do_map(expr): @@ -1390,14 +1322,12 @@ class InverseMassContractor(CSECachingMapperMixin, IdentityMapper): IdentityMapper.map_common_subexpression def map_boundary_pair(self, bp): - from grudge.optemplate import BoundaryPair - return BoundaryPair(self.rec(bp.field), self.rec(bp.bfield), bp.tag) + return sym.BoundaryPair(self.rec(bp.field), self.rec(bp.bfield), bp.tag) def map_operator_binding(self, binding): # we only care about bindings of inverse mass operators - from grudge.optemplate import InverseMassOperator - if isinstance(binding.op, InverseMassOperator): + if isinstance(binding.op, op.InverseMassOperator): iimc = _InnerInverseMassContractor(self) proposed_result = iimc(binding.field) if iimc.extra_operator_count > 1: @@ -1422,11 +1352,9 @@ class ErrorChecker(CSECachingMapperMixin, IdentityMapper): self.mesh = mesh def map_operator_binding(self, expr): - from grudge.optemplate import DiffOperatorBase - - if isinstance(expr.op, DiffOperatorBase): + if isinstance(expr.op, op.DiffOperatorBase): if (self.mesh is not None - and expr.op.xyz_axis >= self.mesh.dimensions): + and expr.op.xyz_axis >= self.mesh.dim): raise ValueError("optemplate tries to differentiate along a " "non-existent axis (e.g. Z in 2D)") @@ -1482,9 +1410,7 @@ class FluxCollector(CSECachingMapperMixin, CollectorMixin, CombineMapper): CombineMapper.map_common_subexpression def map_operator_binding(self, expr): - from grudge.optemplate import FluxOperatorBase - - if isinstance(expr.op, FluxOperatorBase): + if isinstance(expr.op, op.FluxOperatorBase): result = set([expr]) else: result = set() @@ -1545,17 +1471,16 @@ class FluxExchangeCollector(CSECachingMapperMixin, CollectorMixin, CombineMapper class Evaluator(pymbolic.mapper.evaluator.EvaluationMapper): def map_boundary_pair(self, bp): - from grudge.optemplate.primitives import BoundaryPair - return BoundaryPair(self.rec(bp.field), self.rec(bp.bfield), bp.tag) + return sym.BoundaryPair(self.rec(bp.field), self.rec(bp.bfield), bp.tag) # }}} -# {{{ boundary combiner (used by CUDA backend) +# {{{ boundary combiner class BoundaryCombiner(CSECachingMapperMixin, IdentityMapper): """Combines inner fluxes and boundary fluxes into a single, whole-domain operator of type - :class:`grudge.optemplate.operators.WholeDomainFluxOperator`. + :class:`grudge.symbolic.operators.WholeDomainFluxOperator`. """ def __init__(self, mesh): self.mesh = mesh @@ -1564,8 +1489,7 @@ class BoundaryCombiner(CSECachingMapperMixin, IdentityMapper): IdentityMapper.map_common_subexpression def gather_one_wdflux(self, expressions): - from grudge.optemplate.primitives import OperatorBinding, BoundaryPair - from grudge.optemplate.operators import WholeDomainFluxOperator + from grudge.symbolic.primitives import OperatorBinding, BoundaryPair interiors = [] boundaries = [] @@ -1578,9 +1502,8 @@ class BoundaryCombiner(CSECachingMapperMixin, IdentityMapper): rest = [] for ch in expressions: - from grudge.optemplate.operators import FluxOperatorBase if (isinstance(ch, OperatorBinding) - and isinstance(ch.op, FluxOperatorBase)): + and isinstance(ch.op, op.FluxOperatorBase)): skip = False my_is_lift = ch.op.is_lift @@ -1591,10 +1514,7 @@ class BoundaryCombiner(CSECachingMapperMixin, IdentityMapper): if is_lift != my_is_lift: skip = True - from grudge.optemplate.operators import \ - QuadratureFluxOperatorBase - - if isinstance(ch.op, QuadratureFluxOperatorBase): + if isinstance(ch.op, op.QuadratureFluxOperatorBase): my_quad_tag = ch.op.quadrature_tag else: my_quad_tag = None @@ -1612,18 +1532,18 @@ class BoundaryCombiner(CSECachingMapperMixin, IdentityMapper): if isinstance(ch.field, BoundaryPair): bpair = self.rec(ch.field) if self.mesh.tag_to_boundary.get(bpair.tag, []): - boundaries.append(WholeDomainFluxOperator.BoundaryInfo( + boundaries.append(op.WholeDomainFluxOperator.BoundaryInfo( flux_expr=ch.op.flux, bpair=bpair)) else: - interiors.append(WholeDomainFluxOperator.InteriorInfo( + interiors.append(op.WholeDomainFluxOperator.InteriorInfo( flux_expr=ch.op.flux, field_expr=self.rec(ch.field))) else: rest.append(ch) if interiors or boundaries: - wdf = WholeDomainFluxOperator( + wdf = op.WholeDomainFluxOperator( is_lift=is_lift, interiors=interiors, boundaries=boundaries, @@ -1633,8 +1553,7 @@ class BoundaryCombiner(CSECachingMapperMixin, IdentityMapper): return wdf, rest def map_operator_binding(self, expr): - from grudge.optemplate.operators import FluxOperatorBase - if isinstance(expr.op, FluxOperatorBase): + if isinstance(expr.op, op.FluxOperatorBase): wdf, rest = self.gather_one_wdflux([expr]) assert not rest return wdf diff --git a/grudge/symbolic/mappers/bc_to_flux.py b/grudge/symbolic/mappers/bc_to_flux.py index aa2fe090fca46307dad4ff2db2d58634c9dcd927..c3307688c584a1b6dbae79d4e997374f0e99c35c 100644 --- a/grudge/symbolic/mappers/bc_to_flux.py +++ b/grudge/symbolic/mappers/bc_to_flux.py @@ -27,9 +27,11 @@ THE SOFTWARE. from pytools import memoize_method from pymbolic.mapper import CSECachingMapperMixin -from grudge.optemplate.mappers import ( +from grudge.symbolic.mappers import ( IdentityMapper, DependencyMapper, CombineMapper, OperatorReducerMixin) +from grudge import sym +from grudge import sym_flux class ExpensiveBoundaryOperatorDetector(CombineMapper): @@ -41,21 +43,16 @@ class ExpensiveBoundaryOperatorDetector(CombineMapper): return False def map_operator_binding(self, expr): - from grudge.optemplate import (BoundarizeOperator, - FluxExchangeOperator, - QuadratureGridUpsampler, - QuadratureBoundaryGridUpsampler) - - if isinstance(expr.op, BoundarizeOperator): + if isinstance(expr.op, sym.BoundarizeOperator): return False - elif isinstance(expr.op, FluxExchangeOperator): + elif isinstance(expr.op, sym.FluxExchangeOperator): # FIXME: Duplication of these is an even bigger problem! return True elif isinstance(expr.op, ( - QuadratureGridUpsampler, - QuadratureBoundaryGridUpsampler)): + sym.QuadratureGridUpsampler, + sym.QuadratureBoundaryGridUpsampler)): return True else: @@ -98,12 +95,10 @@ class BCToFluxRewriter(CSECachingMapperMixin, IdentityMapper): ExpensiveBoundaryOperatorDetector() def map_operator_binding(self, expr): - from grudge.optemplate.operators import FluxOperatorBase - from grudge.optemplate.primitives import BoundaryPair - from grudge.flux import FluxSubstitutionMapper, FieldComponent + from grudge.symbolic.flux.mappers import FluxSubstitutionMapper - if not (isinstance(expr.op, FluxOperatorBase) - and isinstance(expr.field, BoundaryPair)): + if not (isinstance(expr.op, sym.FluxOperatorBase) + and isinstance(expr.field, sym.BoundaryPair)): return IdentityMapper.map_operator_binding(self, expr) bpair = expr.field @@ -180,7 +175,7 @@ class BCToFluxRewriter(CSECachingMapperMixin, IdentityMapper): self.expensive_bdry_op_detector(expr.child) if has_expensive_operators: - return FieldComponent( + return sym_flux.FieldComponent( self.register_boundary_expr(expr), is_interior=False) else: @@ -190,7 +185,7 @@ class BCToFluxRewriter(CSECachingMapperMixin, IdentityMapper): raise RuntimeError("Your operator template contains a flux normal. " "You may find this confusing, but you can't do that. " "It turns out that you need to use " - "grudge.optemplate.make_normal() for normals in boundary " + "grudge.sym.normal() for normals in boundary " "terms of operator templates.") def map_normal_component(self, expr): @@ -199,33 +194,27 @@ class BCToFluxRewriter(CSECachingMapperMixin, IdentityMapper): "do not agree about boundary tag: %s vs %s" % (expr.boundary_tag, bpair.tag)) - from grudge.flux import Normal - return Normal(expr.axis) + return sym_flux.Normal(expr.axis) def map_variable(self, expr): - return FieldComponent( + return sym_flux.FieldComponent( self.register_boundary_expr(expr), is_interior=False) map_subscript = map_variable def map_operator_binding(self, expr): - from grudge.optemplate import (BoundarizeOperator, - FluxExchangeOperator, - QuadratureGridUpsampler, - QuadratureBoundaryGridUpsampler) - - if isinstance(expr.op, BoundarizeOperator): + if isinstance(expr.op, sym.BoundarizeOperator): if expr.op.tag != bpair.tag: raise RuntimeError("BoundarizeOperator and BoundaryPair " "do not agree about boundary tag: %s vs %s" % (expr.op.tag, bpair.tag)) - return FieldComponent( + return sym_flux.FieldComponent( self.register_volume_expr(expr.field), is_interior=True) - elif isinstance(expr.op, FluxExchangeOperator): + elif isinstance(expr.op, sym.FluxExchangeOperator): from grudge.mesh import TAG_RANK_BOUNDARY op_tag = TAG_RANK_BOUNDARY(expr.op.rank) if bpair.tag != op_tag: @@ -233,53 +222,51 @@ class BCToFluxRewriter(CSECachingMapperMixin, IdentityMapper): "FluxExchangeOperator do not agree about " "boundary tag: %s vs %s" % (op_tag, bpair.tag)) - return FieldComponent( + return sym_flux.FieldComponent( self.register_boundary_expr(expr), is_interior=False) - elif isinstance(expr.op, QuadratureBoundaryGridUpsampler): + elif isinstance(expr.op, sym.QuadratureBoundaryGridUpsampler): if bpair.tag != expr.op.boundary_tag: raise RuntimeError("BoundarizeOperator " "and QuadratureBoundaryGridUpsampler " "do not agree about boundary tag: %s vs %s" % (expr.op.boundary_tag, bpair.tag)) - return FieldComponent( + return sym_flux.FieldComponent( self.register_boundary_expr(expr), is_interior=False) - elif isinstance(expr.op, QuadratureGridUpsampler): + elif isinstance(expr.op, sym.QuadratureGridUpsampler): # We're invoked before operator specialization, so we may # see these instead of QuadratureBoundaryGridUpsampler. - return FieldComponent( + return sym_flux.FieldComponent( self.register_boundary_expr(expr), is_interior=False) else: raise RuntimeError("Found '%s' in a boundary term. " - "To the best of my knowledge, no grudge operator applies " - "directly to boundary data, so this is likely in error." - % expr.op) + "To the best of my knowledge, no grudge operator applies " + "directly to boundary data, so this is likely in error." + % expr.op) def map_flux_exchange(self, expr): - return FieldComponent( + return sym_flux.FieldComponent( self.register_boundary_expr(expr), is_interior=False) # }}} - from grudge.tools import is_obj_array + from pytools.obj_array import is_obj_array if not is_obj_array(vol_field): vol_field = [vol_field] mbfeef = MaxBoundaryFluxEvaluableExpressionFinder(list(vol_field), self.expensive_bdry_op_detector) - #from grudge.optemplate.tools import pretty - #print pretty(bdry_field) - #raw_input("YO") + new_bdry_field = mbfeef(bdry_field) # Step II: Substitute the new_bdry_field into the flux. def sub_bdry_into_flux(expr): - if isinstance(expr, FieldComponent) and not expr.is_interior: + if isinstance(expr, sym_flux.FieldComponent) and not expr.is_interior: if expr.index == 0 and not is_obj_array(bdry_field): return new_bdry_field else: @@ -289,12 +276,13 @@ class BCToFluxRewriter(CSECachingMapperMixin, IdentityMapper): new_flux = FluxSubstitutionMapper(sub_bdry_into_flux)(flux) - from grudge.tools import is_zero, make_obj_array + from grudge.tools import is_zero + from pytools.obj_array import make_obj_array if is_zero(new_flux): return 0 else: return type(expr.op)(new_flux, *expr.op.__getinitargs__()[1:])( - BoundaryPair( + sym.BoundaryPair( make_obj_array([self.rec(e) for e in mbfeef.vol_expr_list]), make_obj_array([self.rec(e) for e in mbfeef.bdry_expr_list]), bpair.tag)) diff --git a/grudge/symbolic/mappers/type_inference.py b/grudge/symbolic/mappers/type_inference.py index a88cfe5a50ec144782189e56f4854525b7e52ebf..0d77d71a1268a0271ac1e7e69923be7766aeac00 100644 --- a/grudge/symbolic/mappers/type_inference.py +++ b/grudge/symbolic/mappers/type_inference.py @@ -25,6 +25,7 @@ THE SOFTWARE. """ import pymbolic.mapper +from grudge import sym # {{{ representation tags @@ -72,7 +73,7 @@ class QuadratureRepresentation(object): class type_info: """These classes represent various bits and pieces of information that - we may deduce about expressions in our optemplate. + we may deduce about expressions in our symbolic operator. """ # serves only as a namespace, thus lower case @@ -94,10 +95,9 @@ class type_info: if u_s_o is NotImplemented: if u_o_s is NotImplemented: if expr is not None: - from grudge.optemplate.tools import pretty raise TypeError("types '%s' and '%s' for '%s' " "cannot be unified" % (self, other, - pretty(expr))) + sym.pretty(expr))) else: raise TypeError("types '%s' and '%s' cannot be unified" % (self, other)) @@ -461,74 +461,49 @@ class TypeInferrer(pymbolic.mapper.RecursiveMapper): children=expr.parameters) def map_operator_binding(self, expr, typedict): - from grudge.optemplate.operators import ( - NodalReductionOperator, - - DiffOperatorBase, - ReferenceDiffOperatorBase, - - MassOperatorBase, - ReferenceMassOperatorBase, - - ElementwiseMaxOperator, - BoundarizeOperator, FluxExchangeOperator, - FluxOperatorBase, - QuadratureGridUpsampler, QuadratureBoundaryGridUpsampler, - QuadratureInteriorFacesGridUpsampler, - - MassOperator, - ReferenceMassOperator, - ReferenceQuadratureMassOperator, - - StiffnessTOperator, - ReferenceStiffnessTOperator, - ReferenceQuadratureStiffnessTOperator, - - ElementwiseLinearOperator) - - if isinstance(expr.op, NodalReductionOperator): + if isinstance(expr.op, sym.NodalReductionOperator): typedict[expr.field] = type_info.KnownVolume() self.rec(expr.field, typedict) return type_info.Scalar() elif isinstance(expr.op, - (ReferenceQuadratureStiffnessTOperator, - ReferenceQuadratureMassOperator)): + (sym.ReferenceQuadratureStiffnessTOperator, + sym.ReferenceQuadratureMassOperator)): typedict[expr.field] = type_info.VolumeVector( QuadratureRepresentation(expr.op.quadrature_tag)) self.rec(expr.field, typedict) return type_info.VolumeVector(NodalRepresentation()) elif isinstance(expr.op, - (ReferenceStiffnessTOperator, StiffnessTOperator)): + (sym.ReferenceStiffnessTOperator, sym.StiffnessTOperator)): # stiffness_T can be specialized for quadrature by OperatorSpecializer typedict[expr.field] = type_info.KnownVolume() self.rec(expr.field, typedict) return type_info.VolumeVector(NodalRepresentation()) elif isinstance(expr.op, - (ReferenceMassOperator, MassOperator)): + (sym.ReferenceMassOperator, sym.MassOperator)): # mass can be specialized for quadrature by OperatorSpecializer typedict[expr.field] = type_info.KnownVolume() self.rec(expr.field, typedict) return type_info.VolumeVector(NodalRepresentation()) elif isinstance(expr.op, ( - DiffOperatorBase, - ReferenceDiffOperatorBase, - MassOperatorBase, - ReferenceMassOperatorBase)): + sym.DiffOperatorBase, + sym.ReferenceDiffOperatorBase, + sym.MassOperatorBase, + sym.ReferenceMassOperatorBase)): # all other operators are purely nodal typedict[expr.field] = type_info.VolumeVector(NodalRepresentation()) self.rec(expr.field, typedict) return type_info.VolumeVector(NodalRepresentation()) - elif isinstance(expr.op, ElementwiseMaxOperator): + elif isinstance(expr.op, sym.ElementwiseMaxOperator): typedict[expr.field] = typedict[expr].unify( type_info.KnownVolume(), expr.field) return self.rec(expr.field, typedict) - elif isinstance(expr.op, BoundarizeOperator): + elif isinstance(expr.op, sym.BoundarizeOperator): # upward propagation: argument has same rep tag as result typedict[expr.field] = type_info.KnownVolume().unify( extract_representation(type_info), expr.field) @@ -539,12 +514,11 @@ class TypeInferrer(pymbolic.mapper.RecursiveMapper): return type_info.KnownBoundary(expr.op.tag).unify( extract_representation(typedict[expr.field]), expr) - elif isinstance(expr.op, FluxExchangeOperator): + elif isinstance(expr.op, sym.FluxExchangeOperator): raise NotImplementedError - elif isinstance(expr.op, FluxOperatorBase): + elif isinstance(expr.op, sym.FluxOperatorBase): from pytools.obj_array import with_object_array_or_scalar - from grudge.optemplate.primitives import BoundaryPair repr_tag_cell = [type_info.no_type] @@ -554,7 +528,7 @@ class TypeInferrer(pymbolic.mapper.RecursiveMapper): repr_tag_cell[0] = extract_representation( self.rec(flux_arg, typedict)) - if isinstance(expr.field, BoundaryPair): + if isinstance(expr.field, sym.BoundaryPair): def process_bdry_flux_arg(flux_arg): typedict[flux_arg] = type_info.KnownBoundary(bpair.tag) \ .unify(repr_tag_cell[0], flux_arg) @@ -570,21 +544,21 @@ class TypeInferrer(pymbolic.mapper.RecursiveMapper): return type_info.VolumeVector(NodalRepresentation()) - elif isinstance(expr.op, QuadratureGridUpsampler): + elif isinstance(expr.op, sym.QuadratureGridUpsampler): typedict[expr.field] = extract_domain(typedict[expr]) self.rec(expr.field, typedict) return type_info.KnownRepresentation( QuadratureRepresentation(expr.op.quadrature_tag))\ .unify(extract_domain(typedict[expr.field]), expr) - elif isinstance(expr.op, QuadratureInteriorFacesGridUpsampler): + elif isinstance(expr.op, sym.QuadratureInteriorFacesGridUpsampler): typedict[expr.field] = type_info.VolumeVector( NodalRepresentation()) self.rec(expr.field, typedict) return type_info.InteriorFacesVector( QuadratureRepresentation(expr.op.quadrature_tag)) - elif isinstance(expr.op, QuadratureBoundaryGridUpsampler): + elif isinstance(expr.op, sym.QuadratureBoundaryGridUpsampler): typedict[expr.field] = type_info.BoundaryVector( expr.op.boundary_tag, NodalRepresentation()) self.rec(expr.field, typedict) @@ -592,7 +566,7 @@ class TypeInferrer(pymbolic.mapper.RecursiveMapper): expr.op.boundary_tag, QuadratureRepresentation(expr.op.quadrature_tag)) - elif isinstance(expr.op, ElementwiseLinearOperator): + elif isinstance(expr.op, sym.ElementwiseLinearOperator): typedict[expr.field] = type_info.VolumeVector(NodalRepresentation()) self.rec(expr.field, typedict) return type_info.VolumeVector(NodalRepresentation()) diff --git a/grudge/symbolic/operators.py b/grudge/symbolic/operators.py index 8b7b635ce6e37f8a33c76744247195405c30c22b..997b171818227fa767f2a355b977752631ce6b4d 100644 --- a/grudge/symbolic/operators.py +++ b/grudge/symbolic/operators.py @@ -778,4 +778,72 @@ class WholeDomainFluxOperator(pymbolic.primitives.AlgebraicLeaf): # }}} +# {{{ convenience functions for operator creation + +def get_flux_operator(flux): + """Return a flux operator that can be multiplied with + a volume field to obtain the interior fluxes + or with a :class:`BoundaryPair` to obtain the lifted boundary + flux. + """ + from pytools.obj_array import is_obj_array + from grudge.symbolic.operators import VectorFluxOperator, FluxOperator + + if is_obj_array(flux): + return VectorFluxOperator(flux) + else: + return FluxOperator(flux) + + +def nabla(dim): + from pytools.obj_array import make_obj_array + return make_obj_array( + [DifferentiationOperator(i) for i in range(dim)]) + + +def minv_stiffness_t(dim): + from pytools.obj_array import make_obj_array + return make_obj_array( + [MInvSTOperator(i) for i in range(dim)]) + + +def stiffness(dim): + from pytools.obj_array import make_obj_array + return make_obj_array( + [StiffnessOperator(i) for i in range(dim)]) + + +def stiffness_t(dim): + from pytools.obj_array import make_obj_array + return make_obj_array( + [StiffnessTOperator(i) for i in range(dim)]) + + +def integral(arg): + import grudge.symbolic as sym + return sym.NodalSum()(sym.MassOperator()(sym.Ones())*arg) + + +def norm(p, arg): + """ + :arg arg: is assumed to be a vector, i.e. have shape ``(n,)``. + """ + import grudge.symbolic as sym + + if p == 2: + comp_norm_squared = sym.NodalSum()( + sym.CFunction("fabs")( + arg * sym.MassOperator()(arg))) + return sym.CFunction("sqrt")(sum(comp_norm_squared)) + + elif p == np.Inf: + comp_norm = sym.NodalMax()(sym.CFunction("fabs")(arg)) + from pymbolic.primitives import Max + return reduce(Max, comp_norm) + + else: + return sum(sym.NodalSum()(sym.CFunction("fabs")(arg)**p))**(1/p) + +# }}} + # vim: foldmethod=marker diff --git a/grudge/symbolic/primitives.py b/grudge/symbolic/primitives.py index 2cfdba75f09e07d92c12f7c8f737fa7ad92521cc..1c8f8a58d22cd7dfc0c4207542610e1cbe90c6f5 100644 --- a/grudge/symbolic/primitives.py +++ b/grudge/symbolic/primitives.py @@ -26,15 +26,13 @@ THE SOFTWARE. from six.moves import range -import numpy +import numpy as np import pymbolic.primitives from meshmode.mesh import BTAG_ALL from pymbolic.primitives import ( # noqa make_common_subexpression, If, Comparison) -from pytools import MovedFunctionDeprecationWrapper - class LeafBase(pymbolic.primitives.AlgebraicLeaf): def stringifier(self): @@ -100,14 +98,14 @@ class OperatorBinding(LeafBase): return self.op, self.field def is_equal(self, other): - from grudge.tools import field_equal + from pytools.obj_array import obj_array_equal return (other.__class__ == self.__class__ and other.op == self.op - and field_equal(other.field, self.field)) + and obj_array_equal(other.field, self.field)) def get_hash(self): - from grudge.tools import hashable_field - return hash((self.__class__, self.op, hashable_field(self.field))) + from pytools.obj_array import obj_array_to_hashable + return hash((self.__class__, self.op, obj_array_to_hashable(self.field))) class PrioritizedSubexpression(pymbolic.primitives.CommonSubexpression): @@ -146,18 +144,18 @@ class BoundaryPair(LeafBase): return (self.field, self.bfield, self.tag) def get_hash(self): - from grudge.tools import hashable_field + from pytools.obj_array import obj_array_to_hashable return hash((self.__class__, - hashable_field(self.field), - hashable_field(self.bfield), + obj_array_to_hashable(self.field), + obj_array_to_hashable(self.bfield), self.tag)) def is_equal(self, other): - from grudge.tools import field_equal + from pytools.obj_array import obj_array_equal return (self.__class__ == other.__class__ - and field_equal(other.field, self.field) - and field_equal(other.bfield, self.bfield) + and obj_array_equal(other.field, self.field) + and obj_array_equal(other.bfield, self.bfield) and other.tag == self.tag) # }}} @@ -187,7 +185,7 @@ class NodeCoordinateComponent(LeafBase): def nodes(dim, quadrature_tag=None): - return numpy.array([NodeCoordinateComponent(i, quadrature_tag) + return np.array([NodeCoordinateComponent(i, quadrature_tag) for i in range(dim)], dtype=object) @@ -204,11 +202,9 @@ class BoundaryNormalComponent(LeafBase): def normal(tag, dimensions): - return numpy.array([BoundaryNormalComponent(tag, i) + return np.array([BoundaryNormalComponent(tag, i) for i in range(dimensions)], dtype=object) -make_normal = MovedFunctionDeprecationWrapper(normal) - class GeometricFactorBase(LeafBase): def __init__(self, quadrature_tag): diff --git a/grudge/symbolic/tools.py b/grudge/symbolic/tools.py index d5df55dd3ae4354c9eb31392e2a4aedc3d39f5fa..843a538e56b6e3a2b67fe28a6e1a681ee5b589ab 100644 --- a/grudge/symbolic/tools.py +++ b/grudge/symbolic/tools.py @@ -25,108 +25,19 @@ THE SOFTWARE. """ -from six.moves import range, reduce +from six.moves import range import numpy as np -import pymbolic.primitives # noqa -from decorator import decorator -# {{{ convenience functions for optemplate creation - -def get_flux_operator(flux): - """Return a flux operator that can be multiplied with - a volume field to obtain the interior fluxes - or with a :class:`BoundaryPair` to obtain the lifted boundary - flux. - """ - from grudge.tools import is_obj_array - from grudge.symbolic import VectorFluxOperator, FluxOperator - - if is_obj_array(flux): - return VectorFluxOperator(flux) - else: - return FluxOperator(flux) - - -def make_nabla(dim): - from grudge.tools import make_obj_array - from grudge.symbolic import DifferentiationOperator - return make_obj_array( - [DifferentiationOperator(i) for i in range(dim)]) - - -def make_minv_stiffness_t(dim): - from grudge.tools import make_obj_array - from grudge.symbolic import MInvSTOperator - return make_obj_array( - [MInvSTOperator(i) for i in range(dim)]) - - -def make_stiffness(dim): - from grudge.tools import make_obj_array - from grudge.symbolic import StiffnessOperator - return make_obj_array( - [StiffnessOperator(i) for i in range(dim)]) - - -def make_stiffness_t(dim): - from grudge.tools import make_obj_array - from grudge.symbolic import StiffnessTOperator - return make_obj_array( - [StiffnessTOperator(i) for i in range(dim)]) - - -def integral(arg): - import grudge.symbolic as sym - return sym.NodalSum()(sym.MassOperator()(sym.Ones())*arg) - - -def norm(p, arg): - """ - :arg arg: is assumed to be a vector, i.e. have shape ``(n,)``. - """ - import grudge.symbolic as sym - - if p == 2: - comp_norm_squared = sym.NodalSum()( - sym.CFunction("fabs")( - arg * sym.MassOperator()(arg))) - return sym.CFunction("sqrt")(sum(comp_norm_squared)) - - elif p == np.Inf: - comp_norm = sym.NodalMax()(sym.CFunction("fabs")(arg)) - from pymbolic.primitives import Max - return reduce(Max, comp_norm) - - else: - return sum(sym.NodalSum()(sym.CFunction("fabs")(arg)**p))**(1/p) - - -def flat_end_sin(x): - from grudge.symbolic.primitives import CFunction - from pymbolic.primitives import IfPositive - from math import pi - return IfPositive(-pi/2-x, - -1, IfPositive(x-pi/2, 1, CFunction("sin")(x))) - - -def smooth_ifpos(crit, right, left, width): - from math import pi - return 0.5*((left+right) - + (right-left)*flat_end_sin( - pi/2/width * crit)) -# }}} - - -# {{{ optemplate tools +# {{{ symbolic operator tools def is_scalar(expr): - from grudge.symbolic import ScalarParameter - return isinstance(expr, (int, float, complex, ScalarParameter)) + from grudge import sym + return isinstance(expr, (int, float, complex, sym.ScalarParameter)) def get_flux_dependencies(flux, field, bdry="all"): - from grudge.flux import FluxDependencyMapper, FieldComponent + from grudge.symbolic.flux.mappers import FluxDependencyMapper, FieldComponent in_fields = list(FluxDependencyMapper( include_calls=False)(flux)) @@ -143,8 +54,8 @@ def get_flux_dependencies(flux, field, bdry="all"): return fld from grudge.tools import is_zero - from grudge.symbolic import BoundaryPair - if isinstance(field, BoundaryPair): + from grudge import sym + if isinstance(field, sym.BoundaryPair): for inf in in_fields: if inf.is_interior: if bdry in ["all", "int"]: @@ -165,7 +76,7 @@ def get_flux_dependencies(flux, field, bdry="all"): yield value -def split_optemplate_for_multirate(state_vector, sym_operator, +def split_sym_operator_for_multirate(state_vector, sym_operator, index_groups): class IndexGroupKillerSubstMap: def __init__(self, kill_set): @@ -246,92 +157,14 @@ def ptwise_dot(logdims1, logdims2, a1, a2): # }}} -# {{{ process_optemplate function - -def process_optemplate(optemplate, post_bind_mapper=None, - dumper=lambda name, optemplate: None, mesh=None, - type_hints={}): - - from grudge.symbolic.mappers import ( - OperatorBinder, CommutativeConstantFoldingMapper, - EmptyFluxKiller, InverseMassContractor, DerivativeJoiner, - ErrorChecker, OperatorSpecializer, GlobalToReferenceMapper) - from grudge.symbolic.mappers.bc_to_flux import BCToFluxRewriter - from grudge.symbolic.mappers.type_inference import TypeInferrer - - dumper("before-bind", optemplate) - optemplate = OperatorBinder()(optemplate) - - ErrorChecker(mesh)(optemplate) - - if post_bind_mapper is not None: - dumper("before-postbind", optemplate) - optemplate = post_bind_mapper(optemplate) - - if mesh is not None: - dumper("before-empty-flux-killer", optemplate) - optemplate = EmptyFluxKiller(mesh)(optemplate) - - dumper("before-cfold", optemplate) - optemplate = CommutativeConstantFoldingMapper()(optemplate) - - dumper("before-bc2flux", optemplate) - optemplate = BCToFluxRewriter()(optemplate) - - # Ordering restriction: - # - # - Must run constant fold before first type inference pass, because zeros, - # while allowed, violate typing constraints (because they can't be assigned - # a unique type), and need to be killed before the type inferrer sees them. - # - # - Must run BC-to-flux before first type inferrer run so that zeros in - # flux arguments can be removed. - - dumper("before-specializer", optemplate) - optemplate = OperatorSpecializer( - TypeInferrer()(optemplate, type_hints) - )(optemplate) - - # Ordering restriction: - # - # - Must run OperatorSpecializer before performing the GlobalToReferenceMapper, - # because otherwise it won't differentiate the type of grids (node or quadrature - # grids) that the operators will apply on. - - assert mesh is not None - dumper("before-global-to-reference", optemplate) - optemplate = GlobalToReferenceMapper(mesh.dimensions)(optemplate) - - # Ordering restriction: - # - # - Must specialize quadrature operators before performing inverse mass - # contraction, because there are no inverse-mass-contracted variants of the - # quadrature operators. - - dumper("before-imass", optemplate) - optemplate = InverseMassContractor()(optemplate) - - dumper("before-cfold-2", optemplate) - optemplate = CommutativeConstantFoldingMapper()(optemplate) - - dumper("before-derivative-join", optemplate) - optemplate = DerivativeJoiner()(optemplate) - - dumper("process-optemplate-finished", optemplate) - - return optemplate - -# }}} - - # {{{ pretty printing -def pretty(optemplate): +def pretty(sym_operator): from grudge.symbolic.mappers import PrettyStringifyMapper stringify_mapper = PrettyStringifyMapper() from pymbolic.mapper.stringifier import PREC_NONE - result = stringify_mapper(optemplate, PREC_NONE) + result = stringify_mapper(sym_operator, PREC_NONE) splitter = "="*75 + "\n" @@ -356,32 +189,4 @@ def pretty(optemplate): # }}} -@decorator -def memoize_method_with_obj_array_args(method, instance, *args): - """This decorator manages to memoize functions that - take object arrays (which are mutable, but are assumed - to never change) as arguments. - """ - dicname = "_memoize_dic_"+method.__name__ - - new_args = [] - for arg in args: - if isinstance(arg, np.ndarray) and arg.dtype == object: - new_args.append(tuple(arg)) - else: - new_args.append(arg) - new_args = tuple(new_args) - - try: - return getattr(instance, dicname)[new_args] - except AttributeError: - result = method(instance, *args) - setattr(instance, dicname, {new_args: result}) - return result - except KeyError: - result = method(instance, *args) - getattr(instance, dicname)[new_args] = result - return result - - # vim: foldmethod=marker diff --git a/grudge/tools.py b/grudge/tools.py new file mode 100644 index 0000000000000000000000000000000000000000..cd1cc49f1828ba82f95c0e805262315e5f8e2be6 --- /dev/null +++ b/grudge/tools.py @@ -0,0 +1,38 @@ +"""Miscellaneous helper facilities.""" + +from __future__ import division + +__copyright__ = "Copyright (C) 2007 Andreas Kloeckner" + +__license__ = """ +Permission is hereby granted, free of charge, to any person obtaining a copy +of this software and associated documentation files (the "Software"), to deal +in the Software without restriction, including without limitation the rights +to use, copy, modify, merge, publish, distribute, sublicense, and/or sell +copies of the Software, and to permit persons to whom the Software is +furnished to do so, subject to the following conditions: + +The above copyright notice and this permission notice shall be included in +all copies or substantial portions of the Software. + +THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR +IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY, +FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL THE +AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER +LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM, +OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN +THE SOFTWARE. +""" + +import numpy as np + + +def is_zero(x): + # DO NOT try to replace this with an attempted '== 0' comparison. + # This will become an elementwise numpy operation and not do what + # you want. + + if np.isscalar(x): + return x == 0 + else: + return False