diff --git a/examples/wave/wave-min.py b/examples/wave/wave-min.py index bf2cc843f19314fd13463e1263499150817cba7c..6e2baa1bcc328af322a8e45cc1d04c2f78613a8b 100644 --- a/examples/wave/wave-min.py +++ b/examples/wave/wave-min.py @@ -35,12 +35,19 @@ def main(write_output=True, order=4): cl_ctx = cl.create_some_context() queue = cl.CommandQueue(cl_ctx) - dims = 2 + dims = 3 from meshmode.mesh.generation import generate_regular_rect_mesh mesh = generate_regular_rect_mesh( a=(-0.5,)*dims, b=(0.5,)*dims, - n=(8,)*dims) + n=(16,)*dims) + + if mesh.dim == 2: + dt = 0.04 + elif mesh.dim == 3: + dt = 0.02 + + print("%d elements" % mesh.nelements) discr = Discretization(cl_ctx, mesh, order=order) @@ -77,17 +84,12 @@ def main(write_output=True, order=4): # print(sym.pretty(op.sym_operator())) bound_op = bind(discr, op.sym_operator()) - print(bound_op) + # print(bound_op) # 1/0 def rhs(t, w): return bound_op(queue, t=t, w=w) - if mesh.dim == 2: - dt = 0.04 - elif mesh.dim == 3: - dt = 0.02 - dt_stepper = set_up_rk4("w", dt, fields, rhs) final_t = 10 diff --git a/grudge/execution.py b/grudge/execution.py index e3646b7361d8b1606a1a234fbe6089a51127f8b7..69e4b7a420cdc55f0a90d1188344f9d1d5ecde17 100644 --- a/grudge/execution.py +++ b/grudge/execution.py @@ -81,11 +81,6 @@ class ExecutionMapper(mappers.Evaluator, discr = self.get_discr(expr.dd) return discr.nodes()[expr.axis].with_queue(self.queue) - 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_grudge_variable(self, expr): return self.context[expr.name] @@ -105,13 +100,16 @@ class ExecutionMapper(mappers.Evaluator, return func(*pars) def map_nodal_sum(self, op, field_expr): - return cl.array.sum(self.rec(field_expr)) + # FIXME: Could allow array scalars + return cl.array.sum(self.rec(field_expr)).get()[()] def map_nodal_max(self, op, field_expr): - return cl.array.max(self.rec(field_expr)) + # FIXME: Could allow array scalars + return cl.array.max(self.rec(field_expr)).get()[()] def map_nodal_min(self, op, field_expr): - return cl.array.min(self.rec(field_expr)) + # FIXME: Could allow array scalars + return cl.array.min(self.rec(field_expr)).get()[()] def map_if(self, expr): bool_crit = self.rec(expr.condition) @@ -333,11 +331,34 @@ class ExecutionMapper(mappers.Evaluator, # {{{ code execution functions - def exec_assign(self, insn): + def map_insn_loopy_kernel(self, insn): + kwargs = {} + kdescr = insn.kernel_descriptor + for name, expr in six.iteritems(kdescr.input_mappings): + kwargs[name] = self.rec(expr) + + discr = self.get_discr(kdescr.governing_dd) + for name in kdescr.scalar_args(): + v = kwargs[name] + if isinstance(v, (int, float)): + kwargs[name] = discr.real_dtype.type(v) + elif isinstance(v, complex): + kwargs[name] = discr.complex_dtype.type(v) + elif isinstance(v, np.number): + pass + else: + raise ValueError("unrecognized scalar type for variable '%s': %s" + % (name, type(v))) + + kwargs["grdg_n"] = discr.nnodes + evt, result_dict = kdescr.loopy_kernel(self.queue, **kwargs) + return list(result_dict.items()), [] + + def map_insn_assign(self, insn): return [(name, self.rec(expr)) for name, expr in zip(insn.names, insn.exprs)], [] - def exec_assign_to_discr_scoped(self, insn): + def map_insn_assign_to_discr_scoped(self, insn): assignments = [] for name, expr in zip(insn.names, insn.exprs): value = self.rec(expr) @@ -346,11 +367,11 @@ class ExecutionMapper(mappers.Evaluator, return assignments, [] - def exec_assign_from_discr_scoped(self, insn): + def map_insn_assign_from_discr_scoped(self, insn): return [(insn.name, self.discr._discr_scoped_subexpr_name_to_value[insn.name])], [] - def exec_diff_batch_assign(self, insn): + def map_insn_diff_batch_assign(self, insn): field = self.rec(insn.field) repr_op = insn.operators[0] # FIXME: There's no real reason why differentiation is special, diff --git a/grudge/symbolic/compiler.py b/grudge/symbolic/compiler.py index fc4ad5f016c522c4c362190dca4ae76aedd5cbd5..ac7765ed919f87ed7b8a17099a302263b94ce347 100644 --- a/grudge/symbolic/compiler.py +++ b/grudge/symbolic/compiler.py @@ -31,6 +31,7 @@ from pytools import Record, memoize_method, memoize from grudge import sym import grudge.symbolic.mappers as mappers from pymbolic.primitives import Variable, Subscript +from six.moves import intern # {{{ instructions @@ -38,6 +39,7 @@ from pymbolic.primitives import Variable, Subscript class Instruction(Record): __slots__ = [] priority = 0 + neglect_for_dofdesc_inference = False def get_assignees(self): raise NotImplementedError("no get_assignees in %s" % self.__class__) @@ -48,9 +50,6 @@ class Instruction(Record): def __str__(self): raise NotImplementedError - def get_execution_method(self, exec_mapper): - raise NotImplementedError - def __hash__(self): return id(self) @@ -69,6 +68,60 @@ def _make_dep_mapper(include_subscripts): include_calls="descend_args") +# {{{ loopy kernel instruction + +class LoopyKernelDescriptor(object): + def __init__(self, loopy_kernel, input_mappings, output_mappings, + fixed_arguments, governing_dd): + self.loopy_kernel = loopy_kernel + self.input_mappings = input_mappings + self.output_mappings = output_mappings + self.fixed_arguments = fixed_arguments + self.governing_dd = governing_dd + + @memoize_method + def scalar_args(self): + import loopy as lp + return [arg.name for arg in self.loopy_kernel.args + if isinstance(arg, lp.ValueArg) + and arg.name != "grdg_n"] + + +class LoopyKernelInstruction(Instruction): + comment = "" + scope_indicator = "" + + def __init__(self, kernel_descriptor): + self.kernel_descriptor = kernel_descriptor + + @memoize_method + def get_assignees(self): + return set( + k + for k in self.kernel_descriptor.output_mappings.keys()) + + @memoize_method + def get_dependencies(self): + from pymbolic.primitives import Variable, Subscript + return set( + v + for v in self.kernel_descriptor.input_mappings.values() + if isinstance(v, (Variable, Subscript))) + + def __str__(self): + knl_str = "\n".join( + "%s = %s" % (insn.assignee, insn.expression) + for insn in self.kernel_descriptor.loopy_kernel.instructions) + + knl_str = knl_str.replace("grdg_", "") + + return "{ /* loopy */\n %s\n}" % knl_str.replace("\n", "\n ") + + mapper_method = "map_insn_loopy_kernel" + +# }}} + + class AssignBase(Instruction): comment = "" scope_indicator = "" @@ -110,7 +163,6 @@ class Assign(AssignBase): exprs describes an expression that is not needed beyond this assignment .. attribute:: priority - .. attribute:: is_scalar_valued """ def __init__(self, names, exprs, **kwargs): @@ -143,19 +195,18 @@ class Assign(AssignBase): return deps - def get_execution_method(self, exec_mapper): - return exec_mapper.exec_assign + mapper_method = intern("map_insn_assign") class ToDiscretizationScopedAssign(Assign): scope_indicator = "(to discr)-" - def get_execution_method(self, exec_mapper): - return exec_mapper.exec_assign_to_discr_scoped + mapper_method = intern("map_insn_assign_to_discr_scoped") class FromDiscretizationScopedAssign(AssignBase): scope_indicator = "(discr)-" + neglect_for_dofdesc_inference = True def __init__(self, name, **kwargs): super(FromDiscretizationScopedAssign, self).__init__(name=name, **kwargs) @@ -173,8 +224,7 @@ class FromDiscretizationScopedAssign(AssignBase): def __str__(self): return "%s <-(from discr)" % self.name - def get_execution_method(self, exec_mapper): - return exec_mapper.exec_assign_from_discr_scoped + mapper_method = intern("map_insn_assign_from_discr_scoped") class DiffBatchAssign(Instruction): @@ -212,57 +262,7 @@ class DiffBatchAssign(Instruction): return "\n".join(lines) - def get_execution_method(self, exec_mapper): - return exec_mapper.exec_diff_batch_assign - - -class FluxExchangeBatchAssign(Instruction): - """ - .. attribute:: names - .. attribute:: indices_and_ranks - .. attribute:: rank_to_index_and_name - .. attribute:: arg_fields - """ - - priority = 1 - - def __init__(self, names, indices_and_ranks, arg_fields): - rank_to_index_and_name = {} - for name, (index, rank) in zip( - names, indices_and_ranks): - rank_to_index_and_name.setdefault(rank, []).append( - (index, name)) - - Instruction.__init__(self, - names=names, - indices_and_ranks=indices_and_ranks, - rank_to_index_and_name=rank_to_index_and_name, - arg_fields=arg_fields) - - def get_assignees(self): - return set(self.names) - - @memoize_method - def get_dependencies(self): - dep_mapper = _make_dep_mapper() - result = set() - for fld in self.arg_fields: - result |= dep_mapper(fld) - return result - - def __str__(self): - lines = [] - - lines.append("{") - for n, (index, rank) in zip(self.names, self.indices_and_ranks): - lines.append(" %s <- receive index %s from rank %d [%s]" % ( - n, index, rank, self.arg_fields)) - lines.append("}") - - return "\n".join(lines) - - def get_execution_method(self, exec_mapper): - return exec_mapper.exec_flux_exchange_batch_assign + mapper_method = intern("map_insn_diff_batch_assign") # }}} @@ -283,7 +283,8 @@ def dot_dataflow_graph(code, max_node_label_length=30, node_names[insn] = node_name node_label = str(insn) - if max_node_label_length is not None: + if (max_node_label_length is not None + and not isinstance(insn, LoopyKernelInstruction)): node_label = node_label[:max_node_label_length] if label_wrap_width is not None: @@ -484,8 +485,8 @@ class Code(object): del context[name] done_insns.add(insn) - assignments, new_futures = \ - insn.get_execution_method(exec_mapper)(insn) + mapper_method = getattr(exec_mapper, insn.mapper_method) + assignments, new_futures = mapper_method(insn) if insn is not None: for target, value in assignments: @@ -519,6 +520,7 @@ class Code(object): # }}} # {{{ static schedule execution + class EvaluateFuture(object): """A fake 'instruction' that represents evaluation of a future.""" def __init__(self, future_id): @@ -549,8 +551,8 @@ class Code(object): assignments, new_futures = future() del future else: - assignments, new_futures = \ - insn.get_execution_method(exec_mapper)(insn) + mapper_method = getattr(exec_mapper, insn.mapper_method) + assignments, new_futures = mapper_method(insn) for target, value in assignments: if pre_assign_check is not None: @@ -578,6 +580,375 @@ class Code(object): # }}} +# {{{ assignment aggregration pass + +def aggregate_assignments(inf_mapper, instructions, result, + max_vectors_in_batch_expr): + from pymbolic.primitives import Variable + + # {{{ aggregation helpers + + def get_complete_origins_set(insn, skip_levels=0): + if skip_levels < 0: + skip_levels = 0 + + result = set() + for dep in insn.get_dependencies(): + if isinstance(dep, Variable): + dep_origin = origins_map.get(dep.name, None) + if dep_origin is not None: + if skip_levels <= 0: + result.add(dep_origin) + result |= get_complete_origins_set( + dep_origin, skip_levels-1) + + return result + + var_assignees_cache = {} + + def get_var_assignees(insn): + try: + return var_assignees_cache[insn] + except KeyError: + result = set(Variable(assignee) + for assignee in insn.get_assignees()) + var_assignees_cache[insn] = result + return result + + def aggregate_two_assignments(ass_1, ass_2): + names = ass_1.names + ass_2.names + + from pymbolic.primitives import Variable + deps = (ass_1.get_dependencies() | ass_2.get_dependencies()) \ + - set(Variable(name) for name in names) + + return Assign( + names=names, exprs=ass_1.exprs + ass_2.exprs, + _dependencies=deps, + priority=max(ass_1.priority, ass_2.priority)) + + # }}} + + # {{{ main aggregation pass + + origins_map = dict( + (assignee, insn) + for insn in instructions + for assignee in insn.get_assignees()) + + from pytools import partition + from grudge.symbolic.primitives import DTAG_SCALAR + + unprocessed_assigns, other_insns = partition( + lambda insn: ( + isinstance(insn, Assign) + and not isinstance(insn, ToDiscretizationScopedAssign) + and not isinstance(insn, FromDiscretizationScopedAssign) + and not any( + inf_mapper.infer_for_name(n).domain_tag == DTAG_SCALAR + for n in insn.names)), + instructions) + + # filter out zero-flop-count assigns--no need to bother with those + processed_assigns, unprocessed_assigns = partition( + lambda ass: ass.flop_count() == 0, + unprocessed_assigns) + + # filter out zero assignments + from grudge.tools import is_zero + + i = 0 + + while i < len(unprocessed_assigns): + my_assign = unprocessed_assigns[i] + if any(is_zero(expr) for expr in my_assign.exprs): + processed_assigns.append(unprocessed_assigns.pop(i)) + else: + i += 1 + + # greedy aggregation + while unprocessed_assigns: + my_assign = unprocessed_assigns.pop() + + my_deps = my_assign.get_dependencies() + my_assignees = get_var_assignees(my_assign) + + agg_candidates = [] + for i, other_assign in enumerate(unprocessed_assigns): + other_deps = other_assign.get_dependencies() + other_assignees = get_var_assignees(other_assign) + + if ((my_deps & other_deps + or my_deps & other_assignees + or other_deps & my_assignees) + and my_assign.priority == other_assign.priority): + agg_candidates.append((i, other_assign)) + + did_work = False + + if agg_candidates: + my_indirect_origins = get_complete_origins_set( + my_assign, skip_levels=1) + + for other_assign_index, other_assign in agg_candidates: + if max_vectors_in_batch_expr is not None: + new_assignee_count = len( + set(my_assign.get_assignees()) + | set(other_assign.get_assignees())) + new_dep_count = len( + my_assign.get_dependencies( + each_vector=True) + | other_assign.get_dependencies( + each_vector=True)) + + if (new_assignee_count + new_dep_count + > max_vectors_in_batch_expr): + continue + + other_indirect_origins = get_complete_origins_set( + other_assign, skip_levels=1) + + if (my_assign not in other_indirect_origins and + other_assign not in my_indirect_origins): + did_work = True + + # aggregate the two assignments + new_assignment = aggregate_two_assignments( + my_assign, other_assign) + del unprocessed_assigns[other_assign_index] + unprocessed_assigns.append(new_assignment) + for assignee in new_assignment.get_assignees(): + origins_map[assignee] = new_assignment + + break + + if not did_work: + processed_assigns.append(my_assign) + + externally_used_names = set( + expr + for insn in processed_assigns + other_insns + for expr in insn.get_dependencies()) + + from pytools.obj_array import is_obj_array + if is_obj_array(result): + externally_used_names |= set(expr for expr in result) + else: + externally_used_names |= set([result]) + + def schedule_and_finalize_assignment(ass): + dep_mapper = _make_dep_mapper(include_subscripts=False) + + names_exprs = list(zip(ass.names, ass.exprs)) + + my_assignees = set(name for name, expr in names_exprs) + names_exprs_deps = [ + (name, expr, + set(dep.name for dep in dep_mapper(expr) if + isinstance(dep, Variable)) & my_assignees) + for name, expr in names_exprs] + + ordered_names_exprs = [] + available_names = set() + + while names_exprs_deps: + schedulable = [] + + i = 0 + while i < len(names_exprs_deps): + name, expr, deps = names_exprs_deps[i] + + unsatisfied_deps = deps - available_names + + if not unsatisfied_deps: + schedulable.append((str(expr), name, expr)) + del names_exprs_deps[i] + else: + i += 1 + + # make sure these come out in a constant order + schedulable.sort() + + if schedulable: + for key, name, expr in schedulable: + ordered_names_exprs.append((name, expr)) + available_names.add(name) + else: + raise RuntimeError("aggregation resulted in an " + "impossible assignment") + + return Assign( + names=[name for name, expr in ordered_names_exprs], + exprs=[expr for name, expr in ordered_names_exprs], + do_not_return=[Variable(name) not in externally_used_names + for name, expr in ordered_names_exprs], + priority=ass.priority) + + return [schedule_and_finalize_assignment(ass) + for ass in processed_assigns] + other_insns + + # }}} + +# }}} + + +# {{{ + +def set_once(d, k, v): + try: + v_prev = d[k] + except KeyError: + d[k] = v + else: + assert v_prev == d[k] + + +class ToLoopyExpressionMapper(mappers.IdentityMapper): + def __init__(self, dd_inference_mapper, output_names, temp_names, iname): + self.dd_inference_mapper = dd_inference_mapper + self.output_names = output_names + self.temp_names = temp_names + self.iname = iname + from pymbolic import var + self.iname_expr = var(iname) + + self.input_mappings = {} + self.output_mappings = {} + self.non_scalar_vars = [] + + def map_name(self, name): + dot_idx = name.find(".") + if dot_idx != -1: + return "grdg_sub_%s_%s" % (name[:dot_idx], name[dot_idx+1:]) + else: + return name + + def map_variable_reference(self, name, expr): + from pymbolic import var + dd = self.dd_inference_mapper(expr) + + mapped_name = self.map_name(name) + if name in self.output_names: + set_once(self.output_mappings, name, expr) + else: + set_once(self.input_mappings, mapped_name, expr) + + from grudge.symbolic.primitives import DTAG_SCALAR + if dd.domain_tag == DTAG_SCALAR or name in self.temp_names: + return var(mapped_name) + else: + self.non_scalar_vars.append(name) + return var(mapped_name)[self.iname_expr] + + def map_variable(self, expr): + return self.map_variable_reference(expr.name, expr) + + def map_grudge_variable(self, expr): + return self.map_variable_reference(expr.name, expr) + + def map_subscript(self, expr): + return self.map_variable_reference(expr.aggregate.name, expr) + + def map_call(self, expr): + if isinstance(expr.function, sym.CFunction): + from pymbolic import var + return var(expr.function.name)( + *[self.rec(par) for par in expr.parameters]) + else: + raise NotImplementedError( + "do not know how to map function '%s' into loopy" + % expr.function) + + def map_ones(self, expr): + return 1 + + def map_node_coordinate_component(self, expr): + mapped_name = "grdg_ncc%d" % expr.axis + set_once(self.input_mappings, mapped_name, expr) + + from pymbolic import var + return var(mapped_name)[self.iname_expr] + + def map_common_subexpression(self, expr): + raise ValueError("not expecting CSEs at this stage in the " + "compilation process") + + +class ToLoopyInstructionMapper(object): + def __init__(self, dd_inference_mapper): + self.dd_inference_mapper = dd_inference_mapper + + def map_insn_assign(self, insn): + from grudge.symbolic.primitives import OperatorBinding + if len(insn.exprs) == 1 and isinstance(insn.exprs[0], OperatorBinding): + return insn + + iname = "grdg_i" + size_name = "grdg_n" + + temp_names = [ + name + for name, dnr in zip(insn.names, insn.do_not_return) + if dnr] + + expr_mapper = ToLoopyExpressionMapper( + self.dd_inference_mapper, insn.names, temp_names, iname) + insns = [] + + import loopy as lp + from pymbolic import var + for name, expr, dnr in zip(insn.names, insn.exprs, insn.do_not_return): + insns.append( + lp.Assignment( + expr_mapper(var(name)), + expr_mapper(expr), + temp_var_type=lp.auto if dnr else None)) + + if not expr_mapper.non_scalar_vars: + return insn + + knl = lp.make_kernel( + "{[%s]: 0 <= %s < %s}" % (iname, iname, size_name), + insns, + default_offset=lp.auto) + + knl = lp.set_options(knl, return_dict=True) + knl = lp.split_iname(knl, iname, 128, outer_tag="g.0", inner_tag="l.0") + + from pytools import single_valued + governing_dd = single_valued( + self.dd_inference_mapper(expr) + for expr in insn.exprs) + + return LoopyKernelInstruction( + LoopyKernelDescriptor( + loopy_kernel=knl, + input_mappings=expr_mapper.input_mappings, + output_mappings=expr_mapper.output_mappings, + fixed_arguments={}, + governing_dd=governing_dd) + ) + + def map_insn_assign_to_discr_scoped(self, insn): + return insn + + def map_insn_assign_from_discr_scoped(self, insn): + return insn + + def map_insn_diff_batch_assign(self, insn): + return insn + + +def rewrite_insn_to_loopy_insns(inf_mapper, insn_list): + insn_mapper = ToLoopyInstructionMapper(inf_mapper) + + return [ + getattr(insn_mapper, insn.mapper_method)(insn) + for insn in insn_list] + +# }}} + + # {{{ compiler class CodeGenerationState(Record): @@ -636,17 +1007,27 @@ class OperatorCompiler(mappers.IdentityMapper): # Finally, walk the expression and build the code. result = super(OperatorCompiler, self).__call__(expr, codegen_state) + eval_code = self.eval_code + del self.eval_code + discr_code = self.discr_code + del self.discr_code + + from grudge.symbolic.dofdesc_inference import DOFDescInferenceMapper + inf_mapper = DOFDescInferenceMapper(discr_code + eval_code) + + eval_code = aggregate_assignments( + inf_mapper, eval_code, result, self.max_vectors_in_batch_expr) + + discr_code = rewrite_insn_to_loopy_insns(inf_mapper, discr_code) + eval_code = rewrite_insn_to_loopy_insns(inf_mapper, eval_code) + from pytools.obj_array import make_obj_array return ( - Code(self.discr_code, + Code(discr_code, make_obj_array( [Variable(name) for name in self.discr_scope_names_copied_to_eval])), - Code( - # FIXME: Enable - #self.aggregate_assignments(self.eval_code, result), - self.eval_code, - result)) + Code(eval_code, result)) # }}} @@ -797,215 +1178,6 @@ class OperatorCompiler(mappers.IdentityMapper): # }}} - # {{{ assignment aggregration pass - - def aggregate_assignments(self, instructions, result): - from pymbolic.primitives import Variable - - # {{{ aggregation helpers - - def get_complete_origins_set(insn, skip_levels=0): - if skip_levels < 0: - skip_levels = 0 - - result = set() - for dep in insn.get_dependencies(): - if isinstance(dep, Variable): - dep_origin = origins_map.get(dep.name, None) - if dep_origin is not None: - if skip_levels <= 0: - result.add(dep_origin) - result |= get_complete_origins_set( - dep_origin, skip_levels-1) - - return result - - var_assignees_cache = {} - - def get_var_assignees(insn): - try: - return var_assignees_cache[insn] - except KeyError: - result = set(Variable(assignee) - for assignee in insn.get_assignees()) - var_assignees_cache[insn] = result - return result - - def aggregate_two_assignments(ass_1, ass_2): - names = ass_1.names + ass_2.names - - from pymbolic.primitives import Variable - deps = (ass_1.get_dependencies() | ass_2.get_dependencies()) \ - - set(Variable(name) for name in names) - - return Assign( - names=names, exprs=ass_1.exprs + ass_2.exprs, - _dependencies=deps, - priority=max(ass_1.priority, ass_2.priority)) - - # }}} - - # {{{ main aggregation pass - - origins_map = dict( - (assignee, insn) - for insn in instructions - for assignee in insn.get_assignees()) - - from pytools import partition - unprocessed_assigns, other_insns = partition( - # FIXME: Re-add check for scalar result, exclude - lambda insn: isinstance(insn, Assign), - instructions) - - # filter out zero-flop-count assigns--no need to bother with those - processed_assigns, unprocessed_assigns = partition( - lambda ass: ass.flop_count() == 0, - unprocessed_assigns) - - # filter out zero assignments - from pytools import any - from grudge.tools import is_zero - - i = 0 - - while i < len(unprocessed_assigns): - my_assign = unprocessed_assigns[i] - if any(is_zero(expr) for expr in my_assign.exprs): - processed_assigns.append(unprocessed_assigns.pop()) - else: - i += 1 - - # greedy aggregation - while unprocessed_assigns: - my_assign = unprocessed_assigns.pop() - - my_deps = my_assign.get_dependencies() - my_assignees = get_var_assignees(my_assign) - - agg_candidates = [] - for i, other_assign in enumerate(unprocessed_assigns): - other_deps = other_assign.get_dependencies() - other_assignees = get_var_assignees(other_assign) - - if ((my_deps & other_deps - or my_deps & other_assignees - or other_deps & my_assignees) - and my_assign.priority == other_assign.priority): - agg_candidates.append((i, other_assign)) - - did_work = False - - if agg_candidates: - my_indirect_origins = get_complete_origins_set( - my_assign, skip_levels=1) - - for other_assign_index, other_assign in agg_candidates: - if self.max_vectors_in_batch_expr is not None: - new_assignee_count = len( - set(my_assign.get_assignees()) - | set(other_assign.get_assignees())) - new_dep_count = len( - my_assign.get_dependencies( - each_vector=True) - | other_assign.get_dependencies( - each_vector=True)) - - if (new_assignee_count + new_dep_count - > self.max_vectors_in_batch_expr): - continue - - other_indirect_origins = get_complete_origins_set( - other_assign, skip_levels=1) - - if (my_assign not in other_indirect_origins and - other_assign not in my_indirect_origins): - did_work = True - - # aggregate the two assignments - new_assignment = aggregate_two_assignments( - my_assign, other_assign) - del unprocessed_assigns[other_assign_index] - unprocessed_assigns.append(new_assignment) - for assignee in new_assignment.get_assignees(): - origins_map[assignee] = new_assignment - - break - - if not did_work: - processed_assigns.append(my_assign) - - externally_used_names = set( - expr - for insn in processed_assigns + other_insns - for expr in insn.get_dependencies()) - - from pytools.obj_array import is_obj_array - if is_obj_array(result): - externally_used_names |= set(expr for expr in result) - else: - externally_used_names |= set([result]) - - def schedule_and_finalize_assignment(ass): - dep_mapper = _make_dep_mapper(include_subscripts=False) - - names_exprs = list(zip(ass.names, ass.exprs)) - - my_assignees = set(name for name, expr in names_exprs) - names_exprs_deps = [ - (name, expr, - set(dep.name for dep in dep_mapper(expr) if - isinstance(dep, Variable)) & my_assignees) - for name, expr in names_exprs] - - ordered_names_exprs = [] - available_names = set() - - while names_exprs_deps: - schedulable = [] - - i = 0 - while i < len(names_exprs_deps): - name, expr, deps = names_exprs_deps[i] - - unsatisfied_deps = deps - available_names - - if not unsatisfied_deps: - schedulable.append((str(expr), name, expr)) - del names_exprs_deps[i] - else: - i += 1 - - # make sure these come out in a constant order - schedulable.sort() - - if schedulable: - for key, name, expr in schedulable: - ordered_names_exprs.append((name, expr)) - available_names.add(name) - else: - raise RuntimeError("aggregation resulted in an " - "impossible assignment") - - return self.finalize_multi_assign( - names=[name for name, expr in ordered_names_exprs], - exprs=[expr for name, expr in ordered_names_exprs], - do_not_return=[Variable(name) not in externally_used_names - for name, expr in ordered_names_exprs], - priority=ass.priority) - - return [schedule_and_finalize_assignment(ass) - for ass in processed_assigns] + other_insns - - # }}} - - # }}} - - def finalize_multi_assign(self, names, exprs, do_not_return, priority): - return Assign(names=names, exprs=exprs, - do_not_return=do_not_return, - priority=priority) - # }}} diff --git a/grudge/symbolic/dofdesc_inference.py b/grudge/symbolic/dofdesc_inference.py new file mode 100644 index 0000000000000000000000000000000000000000..9cb543578266fb5fed990f67ae3758c9af3c18c2 --- /dev/null +++ b/grudge/symbolic/dofdesc_inference.py @@ -0,0 +1,224 @@ +from __future__ import division, absolute_import + +__copyright__ = "Copyright (C) 2017 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. +""" + + +# This is purely leaves-to-roots. No need to propagate information in the +# opposite direction. + + +from pymbolic.mapper import RecursiveMapper, CSECachingMapperMixin +from grudge.symbolic.primitives import DOFDesc, DTAG_SCALAR + + +def unify_dofdescs(dd_a, dd_b, expr=None): + if dd_a is None: + assert dd_b is not None + return dd_b + + if expr is not None: + loc_str = "in expression %s" % str(expr) + else: + loc_str = "" + + from grudge.symbolic.primitives import DTAG_SCALAR + if dd_a.domain_tag != dd_b.domain_tag: + if dd_a.domain_tag == DTAG_SCALAR: + return dd_b + elif dd_b.domain_tag == DTAG_SCALAR: + return dd_a + else: + raise ValueError("mismatched domain tags" + loc_str) + + # domain tags match + if dd_a.quadrature_tag != dd_b.quadrature_tag: + raise ValueError("mismatched quadrature tags" + loc_str) + + return dd_a + + +class InferrableMultiAssignment(object): + """An assignemnt 'instruction' which may be used as part of type + inference. + + .. method:: get_assignees(rec) + + :returns: a :class:`set` of names which are assigned values by + this assignment. + + .. method:: infer_dofdescs(rec) + + :returns: a list of ``(name, :class:`grudge.symbolic.primitives.DOFDesc`)`` + tuples, each indicating the value type of the value with *name*. + """ + + # (not a base class--only documents the interface) + + +class DOFDescInferenceMapper(RecursiveMapper, CSECachingMapperMixin): + def __init__(self, assignments, name_to_dofdesc=None, check=True): + """ + :arg assignments: a list of objects adhering to + :class:`InferrableMultiAssignment`. + :returns: an instance of :class:`DOFDescInferenceMapper` + """ + + self.check = check + + self.name_to_assignment = dict( + (name, a) + for a in assignments + if not a.neglect_for_dofdesc_inference + for name in a.get_assignees()) + + if name_to_dofdesc is None: + name_to_dofdesc = {} + else: + name_to_dofdesc = name_to_dofdesc.copy() + + self.name_to_dofdesc = name_to_dofdesc + + def infer_for_name(self, name): + try: + return self.name_to_dofdesc[name] + except KeyError: + a = self.name_to_assignment[name] + + inf_method = getattr(self, a.mapper_method) + for r_name, r_dofdesc in inf_method(a): + assert r_name not in self.name_to_dofdesc + self.name_to_dofdesc[r_name] = r_dofdesc + + return self.name_to_dofdesc[name] + + # {{{ expression mappings + + def map_constant(self, expr): + return DOFDesc(DTAG_SCALAR) + + def map_grudge_variable(self, expr): + return expr.dd + + def map_variable(self, expr): + return self.infer_for_name(expr.name) + + def map_subscript(self, expr): + # FIXME: Subscript has same type as aggregate--a bit weird + return self.rec(expr.aggregate) + + def map_multi_child(self, expr, children): + dofdesc = None + + for ch in children: + dofdesc = unify_dofdescs(dofdesc, self.rec(ch), expr) + + if dofdesc is None: + raise ValueError("no DOFDesc found for expression %s" % expr) + else: + return dofdesc + + def map_sum(self, expr): + return self.map_multi_child(expr, expr.children) + + map_product = map_sum + + def map_quotient(self, expr): + return self.map_multi_child(expr, (expr.numerator, expr.denominator)) + + def map_power(self, expr): + return self.map_multi_child(expr, (expr.base, expr.exponent)) + + def map_if(self, expr): + return self.map_multi_child(expr, [expr.condition, expr.then, expr.else_]) + + def map_comparison(self, expr): + return self.map_multi_child(expr, [expr.left, expr.right]) + + def map_nodal_sum(self, expr, enclosing_prec): + return DOFDesc(DTAG_SCALAR) + + map_nodal_max = map_nodal_sum + map_nodal_min = map_nodal_sum + + def map_operator_binding(self, expr): + operator = expr.op + + if self.check: + op_dd = self.rec(expr.field) + if op_dd != operator.dd_in: + raise ValueError("mismatched input to %s " + "(got: %s, expected: %s)" + " in '%s'" + % ( + type(expr).__name__, + op_dd, expr.dd_in, + str(expr))) + + return operator.dd_out + + def map_ones(self, expr): + return expr.dd + + map_node_coordinate_component = map_ones + + def map_call(self, expr): + arg_dds = [ + self.rec(par) + for par in expr.parameters] + + assert arg_dds + + # FIXME + return arg_dds[0] + + # }}} + + # {{{ instruction mappings + + def map_insn_assign(self, insn): + return [ + (name, self.rec(expr)) + for name, expr in zip(insn.names, insn.exprs) + ] + + map_insn_assign_to_discr_scoped = map_insn_assign + + def map_insn_diff_batch_assign(self, insn): + if self.check: + repr_op = insn.operators[0] + input_dd = self.rec(insn.field) + if input_dd != repr_op.dd_in: + raise ValueError("mismatched input to %s " + "(got: %s, expected: %s)" + % ( + type(insn).__name__, + input_dd, repr_op.dd_in, + )) + + return [ + (name, op.dd_out) + for name, op in zip(insn.names, insn.operators)] + + # }}} + +# vim: foldmethod=marker diff --git a/grudge/symbolic/mappers/__init__.py b/grudge/symbolic/mappers/__init__.py index 917fa61951b5ad48cd5ab4f3288e1fd662b8bd3f..60b489cecb50d6f69c562ecd7a14e4a0e786ad07 100644 --- a/grudge/symbolic/mappers/__init__.py +++ b/grudge/symbolic/mappers/__init__.py @@ -638,6 +638,7 @@ class StringifyMapper(pymbolic.mapper.stringifier.StringifyMapper): # }}} # {{{ reference differentiation + def map_ref_diff(self, expr, enclosing_prec): return "Diffr%d%s" % (expr.rst_axis, self._format_op_dd(expr)) diff --git a/grudge/symbolic/mappers/bc_to_flux.py b/grudge/symbolic/mappers/bc_to_flux.py deleted file mode 100644 index 545743698a673dd02e12c52343fdb9bfa0778837..0000000000000000000000000000000000000000 --- a/grudge/symbolic/mappers/bc_to_flux.py +++ /dev/null @@ -1,288 +0,0 @@ -"""Operator template mapper: BC-to-flux rewriting.""" - -from __future__ import division - -__copyright__ = "Copyright (C) 2008 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 pytools import memoize_method -from pymbolic.mapper import CSECachingMapperMixin -from grudge.symbolic.mappers import ( - IdentityMapper, DependencyMapper, CombineMapper, - OperatorReducerMixin) -from grudge import sym -from grudge import sym_flux - - -class ExpensiveBoundaryOperatorDetector(CombineMapper): - def combine(self, values): - for val in values: - if val: - return True - - return False - - def map_operator_binding(self, expr): - if isinstance(expr.op, sym.RestrictToBoundary): - return False - - elif isinstance(expr.op, sym.FluxExchangeOperator): - # FIXME: Duplication of these is an even bigger problem! - return True - - elif isinstance(expr.op, ( - sym.QuadratureGridUpsampler, - sym.QuadratureInteriorFacesGridUpsampler)): - return True - - 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) - - def map_common_subexpression(self, expr): - # If there are any expensive operators below here, this - # CSE will catch them, so we can easily flux-CSE down to - # here. - - return False - - def map_normal_component(self, expr): - return False - - map_variable = map_normal_component - map_constant = map_normal_component - - @memoize_method - def __call__(self, expr): - return CombineMapper.__call__(self, expr) - - -class BCToFluxRewriter(CSECachingMapperMixin, IdentityMapper): - """Operates on :class:`FluxOperator` instances bound to - :class:`BoundaryPair`. If the boundary pair's *bfield* is an expression of - what's available in the *field*, we can avoid fetching the data for the - explicit boundary condition and just substitute the *bfield* expression - into the flux. This mapper does exactly that. - """ - - map_common_subexpression_uncached = \ - IdentityMapper.map_common_subexpression - - def __init__(self): - self.expensive_bdry_op_detector = \ - ExpensiveBoundaryOperatorDetector() - - def map_operator_binding(self, expr): - from grudge.symbolic.flux.mappers import FluxSubstitutionMapper - - if not (isinstance(expr.op, sym.FluxOperatorBase) - and isinstance(expr.field, sym.BoundaryPair)): - return IdentityMapper.map_operator_binding(self, expr) - - bpair = expr.field - vol_field = bpair.field - bdry_field = bpair.bfield - flux = expr.op.flux - - bdry_dependencies = DependencyMapper( - include_calls="descend_args", - include_operator_bindings=True)(bdry_field) - - vol_dependencies = DependencyMapper( - include_operator_bindings=True)(vol_field) - - vol_bdry_intersection = bdry_dependencies & vol_dependencies - if vol_bdry_intersection: - raise RuntimeError("Variables are being used as both " - "boundary and volume quantities: %s" - % ", ".join(str(v) for v in vol_bdry_intersection)) - - # Step 1: Find maximal flux-evaluable subexpression of boundary field - # in given BoundaryPair. - - class MaxBoundaryFluxEvaluableExpressionFinder( - IdentityMapper, OperatorReducerMixin): - - def __init__(self, vol_expr_list, expensive_bdry_op_detector): - self.vol_expr_list = vol_expr_list - self.vol_expr_to_idx = dict((vol_expr, idx) - for idx, vol_expr in enumerate(vol_expr_list)) - - self.bdry_expr_list = [] - self.bdry_expr_to_idx = {} - - self.expensive_bdry_op_detector = expensive_bdry_op_detector - - # {{{ expression registration - def register_boundary_expr(self, expr): - try: - return self.bdry_expr_to_idx[expr] - except KeyError: - idx = len(self.bdry_expr_to_idx) - self.bdry_expr_to_idx[expr] = idx - self.bdry_expr_list.append(expr) - return idx - - def register_volume_expr(self, expr): - try: - return self.vol_expr_to_idx[expr] - except KeyError: - idx = len(self.vol_expr_to_idx) - self.vol_expr_to_idx[expr] = idx - self.vol_expr_list.append(expr) - return idx - - # }}} - - # {{{ map_xxx routines - - @memoize_method - def map_common_subexpression(self, expr): - # Here we need to decide whether this CSE should be turned into - # a flux CSE or not. This is a good idea if the transformed - # expression only contains "bare" volume or boundary - # expressions. However, as soon as an operator is applied - # somewhere in the subexpression, the CSE should not be touched - # in order to avoid redundant evaluation of that operator. - # - # Observe that at the time of this writing (Feb 2010), the only - # operators that may occur in boundary expressions are - # quadrature-related. - - has_expensive_operators = \ - self.expensive_bdry_op_detector(expr.child) - - if has_expensive_operators: - return sym_flux.FieldComponent( - self.register_boundary_expr(expr), - is_interior=False) - else: - return IdentityMapper.map_common_subexpression(self, expr) - - def map_normal(self, expr): - 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.sym.normal() for normals in boundary " - "terms of operator templates.") - - def map_normal_component(self, expr): - if expr.boundary_tag != bpair.tag: - raise RuntimeError("BoundaryNormalComponent and BoundaryPair " - "do not agree about boundary tag: %s vs %s" - % (expr.boundary_tag, bpair.tag)) - - return sym_flux.Normal(expr.axis) - - def map_variable(self, expr): - return sym_flux.FieldComponent( - self.register_boundary_expr(expr), - is_interior=False) - - map_subscript = map_variable - - def map_operator_binding(self, expr): - if isinstance(expr.op, sym.RestrictToBoundary): - if expr.op.tag != bpair.tag: - raise RuntimeError("RestrictToBoundary and BoundaryPair " - "do not agree about boundary tag: %s vs %s" - % (expr.op.tag, bpair.tag)) - - return sym_flux.FieldComponent( - self.register_volume_expr(expr.field), - is_interior=True) - - 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: - raise RuntimeError("RestrictToBoundary and " - "FluxExchangeOperator do not agree about " - "boundary tag: %s vs %s" - % (op_tag, bpair.tag)) - return sym_flux.FieldComponent( - self.register_boundary_expr(expr), - is_interior=False) - - elif isinstance(expr.op, sym.QuadratureBoundaryGridUpsampler): - if bpair.tag != expr.op.boundary_tag: - raise RuntimeError("RestrictToBoundary " - "and QuadratureBoundaryGridUpsampler " - "do not agree about boundary tag: %s vs %s" - % (expr.op.boundary_tag, bpair.tag)) - return sym_flux.FieldComponent( - self.register_boundary_expr(expr), - is_interior=False) - - elif isinstance(expr.op, sym.QuadratureGridUpsampler): - # We're invoked before operator specialization, so we may - # see these instead of QuadratureBoundaryGridUpsampler. - 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) - - def map_flux_exchange(self, expr): - return sym_flux.FieldComponent( - self.register_boundary_expr(expr), - is_interior=False) - # }}} - - 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) - - 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, sym_flux.FieldComponent) and not expr.is_interior: - if expr.index == 0 and not is_obj_array(bdry_field): - return new_bdry_field - else: - return new_bdry_field[expr.index] - else: - return None - - new_flux = FluxSubstitutionMapper(sub_bdry_into_flux)(flux) - - 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:])( - 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 deleted file mode 100644 index 68b7d9b30532d6c017bf18922c47806d234c9064..0000000000000000000000000000000000000000 --- a/grudge/symbolic/mappers/type_inference.py +++ /dev/null @@ -1,708 +0,0 @@ -"""Operator template type inference.""" - -from __future__ import division - -__copyright__ = "Copyright (C) 2008 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 six -import pymbolic.mapper -from grudge import sym - - -# {{{ representation tags - -class NodalRepresentation(object): - """A tag representing nodal representation. - - Volume and boundary vectors below are represented either nodally or on a - quadrature grid. This tag expresses one of the two. - """ - def __repr__(self): - return "Nodal" - - def __eq__(self, other): - return type(self) == type(other) - - def __ne__(self, other): - return not self.__eq__(other) - - -class QuadratureRepresentation(object): - """A tag representing representation on a quadrature grid tagged with - *quadrature_tag". - - Volume and boundary vectors below are represented either nodally or on a - quadrature grid. This tag expresses one of the two. - """ - def __init__(self, quadrature_tag): - self.quadrature_tag = quadrature_tag - - def __eq__(self, other): - return (type(self) == type(other) - and self.quadrature_tag == other.quadrature_tag) - - def __ne__(self, other): - return not self.__eq__(other) - - def __repr__(self): - return "Quadrature(%r)" % self.quadrature_tag - -# }}} - - -# {{{ type information -------------------------------------------------------- - -class type_info: # noqa - """These classes represent various bits and pieces of information that - we may deduce about expressions in our symbolic operator. - """ - - # serves only as a namespace, thus lower case - - # {{{ generic type info base classes - class TypeInfo(object): - def unify(self, other, expr=None): - """Return a type that can represent both *self* and *other*. - If impossible, raise :exc:`TypeError`. Subtypes should override - :meth:`unify_inner`. - """ - # shortcut - if self == other: - return self - - u_s_o = self.unify_inner(other) - u_o_s = other.unify_inner(self) - - if u_s_o is NotImplemented: - if u_o_s is NotImplemented: - if expr is not None: - raise TypeError("types '%s' and '%s' for '%s' " - "cannot be unified" % (self, other, - sym.pretty(expr))) - else: - raise TypeError("types '%s' and '%s' cannot be unified" - % (self, other)) - else: - return u_o_s - elif u_o_s is NotImplemented: - return u_s_o - - if u_s_o != u_o_s: - raise RuntimeError("types '%s' and '%s' don't agree about " - "their unifier" % (self, other)) - return u_s_o - - def unify_inner(self, other): - """Actual implementation that tries to unify self and other. - May return *NotImplemented* to indicate that the reverse unification - should be tried. This methods is overriden by derived classes. - Derived classes should delegate to base classes if they don't know the - answer. - """ - return NotImplemented - - def __eq__(self, other): - return (type(self) == type(other) - and self.__getinitargs__() == other.__getinitargs__()) - - def __ne__(self, other): - return not self.__eq__(other) - - class StatelessTypeInfo(TypeInfo): - def __getinitargs__(self): - return () - - class FinalType(TypeInfo): - """If a :class:`TypeInfo` instance is also an instance of this class, - no more information can be added about this type. As a result, this - type only unifies with equal instances. - """ - # }}} - - # {{{ simple types: no type, scalar - class NoType(StatelessTypeInfo): - """Represents "nothing known about type".""" - def unify_inner(self, other): - return other - - def __repr__(self): - return "NoType" - - # this singleton should be the only instance ever created of NoType - no_type = NoType() - - class Scalar(StatelessTypeInfo, FinalType): - def __repr__(self): - return "Scalar" - # }}} - - # {{{ tagged type base classes: representation, domain - class VectorRepresentationBase(object): - def __init__(self, repr_tag): - self.repr_tag = repr_tag - - def __getinitargs__(self): - return (self.repr_tag,) - - class VolumeVectorBase(object): - def __getinitargs__(self): - return () - - class InteriorFacesVectorBase(object): - def __getinitargs__(self): - return () - - class BoundaryVectorBase(object): - def __init__(self, boundary_tag): - self.boundary_tag = boundary_tag - - def __getinitargs__(self): - return (self.boundary_tag,) - # }}} - - # {{{ single-aspect-known unification helper types - class KnownVolume(TypeInfo, VolumeVectorBase): - """Type information indicating that this must be a volume vector - of unknown representation. - """ - - def __repr__(self): - return "KnownAsVolume" - - def unify_inner(self, other): - # Unification with KnownRepresentation is handled in KnownRepresentation. - # Here, we only need to unify with VolumeVector. - - if isinstance(other, type_info.VolumeVector): - return other - else: - return type_info.TypeInfo.unify_inner(self, other) - - class KnownInteriorFaces(TypeInfo, InteriorFacesVectorBase): - """Type information indicating that this must be a vector - of interior face values. - - .. note:: - - These vectors are only necessary in a quadrature setting. - """ - - def __repr__(self): - return "KnownAsIntFace" - - def unify_inner(self, other): - # Unification with KnownRepresentation is handled in KnownRepresentation. - # Here, we only need to unify with InteriorFacesVector. - - if isinstance(other, type_info.InteriorFacesVector): - return other - elif isinstance(other, type_info.KnownVolume): - return type_info.VolumeVector(NodalRepresentation()) - elif other == type_info.VolumeVector(NodalRepresentation()): - return other - else: - return type_info.TypeInfo.unify_inner(self, other) - - class KnownBoundary(TypeInfo, BoundaryVectorBase): - """Type information indicating that this must be a boundary vector.""" - - def __repr__(self): - return "KnownAsBoundary(%s)" % self.boundary_tag - - def unify_inner(self, other): - # Unification with KnownRepresentation is handled in KnownRepresentation. - # Here, we only need to unify with BoundaryVector. - - if (isinstance(other, type_info.BoundaryVector) - and self.boundary_tag == other.boundary_tag): - return other - else: - return type_info.TypeInfo.unify_inner(self, other) - - class KnownRepresentation(TypeInfo, VectorRepresentationBase): - """Type information indicating that the representation (see - representation tags, above) is known, but nothing else (e.g. whether - this is a boundary or volume vector). - """ - def __repr__(self): - return "KnownRepresentation(%s)" % self.repr_tag - - def unify_inner(self, other): - if (isinstance(other, type_info.VolumeVector) - and self.repr_tag == other.repr_tag): - return other - elif (isinstance(other, type_info.BoundaryVector) - and self.repr_tag == other.repr_tag): - return other - elif isinstance(other, type_info.KnownVolume): - return type_info.VolumeVector(self.repr_tag) - elif isinstance(other, type_info.KnownInteriorFaces): - if isinstance(self.repr_tag, NodalRepresentation): - return type_info.VolumeVector(self.repr_tag) - else: - return type_info.InteriorFacesVector(self.repr_tag) - elif isinstance(other, type_info.KnownBoundary): - return type_info.BoundaryVector(other.boundary_tag, self.repr_tag) - else: - return type_info.TypeInfo.unify_inner(self, other) - # }}} - - # {{{ fully specified grudge data types - class VolumeVector(FinalType, VectorRepresentationBase, VolumeVectorBase): - def __repr__(self): - return "Volume(%s)" % self.repr_tag - - class InteriorFacesVector(FinalType, VectorRepresentationBase, - InteriorFacesVectorBase): - def __init__(self, repr_tag): - if not isinstance(repr_tag, QuadratureRepresentation): - raise TypeError("InteriorFacesVector is not usable with non-" - "quadrature representations") - type_info.VectorRepresentationBase.__init__(self, repr_tag) - - def __repr__(self): - return "InteriorFaces(%s)" % self.repr_tag - - class BoundaryVector(FinalType, BoundaryVectorBase, - VectorRepresentationBase): - def __init__(self, boundary_tag, repr_tag): - type_info.BoundaryVectorBase.__init__(self, boundary_tag) - type_info.VectorRepresentationBase.__init__(self, repr_tag) - - def __repr__(self): - return "Boundary(%s, %s)" % (self.boundary_tag, self.repr_tag) - - def __getinitargs__(self): - return (self.boundary_tag, self.repr_tag) - # }}} - - -# {{{ aspect extraction functions - -def extract_representation(ti): - try: - own_repr_tag = ti.repr_tag - except AttributeError: - return type_info.no_type - else: - return type_info.KnownRepresentation(own_repr_tag) - - -def extract_domain(ti): - if isinstance(ti, type_info.VolumeVectorBase): - return type_info.KnownVolume() - elif isinstance(ti, type_info.BoundaryVectorBase): - return type_info.KnownBoundary(ti.boundary_tag) - else: - return type_info.no_type - -# }}} - -# }}} - - -# {{{ TypeDict helper type - -class TypeDict(object): - def __init__(self, hints): - self.container = hints.copy() - self.change_flag = False - - def __getitem__(self, expr): - try: - return self.container[expr] - except KeyError: - return type_info.no_type - - def __setitem__(self, expr, new_tp): - if new_tp is type_info.no_type: - return - - try: - old_tp = self.container[expr] - except KeyError: - self.container[expr] = new_tp - self.change_flag = True - else: - tp = old_tp.unify(new_tp, expr) - if tp != old_tp: - self.change_flag = True - self.container[expr] = tp - - def items(self): - return six.iteritems(self.container) - - iteritems = items - -# }}} - - -# {{{ type inference mapper - -class TypeInferrer(pymbolic.mapper.RecursiveMapper): - def __init__(self): - self.cse_last_results = {} - - def __call__(self, expr, type_hints={}): - typedict = TypeDict(type_hints) - - while True: - typedict.change_flag = False - - def infer_for_expr(expr): - tp = pymbolic.mapper.RecursiveMapper.__call__(self, expr, typedict) - typedict[expr] = tp - - # Numpy arrays occur either at the top level or in flux - # expressions. This code handles the top level case. - from pytools.obj_array import with_object_array_or_scalar - with_object_array_or_scalar(infer_for_expr, expr) - - if not typedict.change_flag: - # nothing has changed any more, type information has 'converged' - break - - # check that type inference completed successfully - for expr, tp in six.iteritems(typedict): - if not isinstance(tp, type_info.FinalType): - raise RuntimeError("type inference was unable to deduce " - "complete type information for '%s' (only '%s')" - % (expr, tp)) - - return typedict - - def rec(self, expr, typedict): - tp = pymbolic.mapper.RecursiveMapper.rec(self, expr, typedict) - typedict[expr] = tp - return tp - - # Information needs to propagate upward (toward the leaves) *and* - # downward (toward the roots) in the expression tree. - - # {{{ base cases - def infer_for_children(self, expr, typedict, children): - # This routine allows scalars among children and treats them as - # not type-changing - - tp = typedict[expr] - - non_scalar_exprs = [] - - for child in children: - if tp is type_info.no_type: - tp = self.rec(child, typedict) - if isinstance(tp, type_info.Scalar): - tp = type_info.no_type - else: - non_scalar_exprs.append(child) - else: - other_tp = self.rec(child, typedict) - - if not isinstance(other_tp, type_info.Scalar): - non_scalar_exprs.append(child) - tp = tp.unify(other_tp, child) - - for child in non_scalar_exprs: - typedict[child] = tp - - if not non_scalar_exprs: - tp = type_info.Scalar() - - return tp - - # }}} - - def map_sum(self, expr, typedict): - return self.infer_for_children(expr, typedict, expr.children) - - def map_product(self, expr, typedict): - return self.infer_for_children(expr, typedict, expr.children) - - def map_quotient(self, expr, typedict): - return self.infer_for_children(expr, typedict, - children=[expr.numerator, expr.denominator]) - - def map_power(self, expr, typedict): - return self.infer_for_children(expr, typedict, - children=[expr.base, expr.exponent]) - - def map_if(self, expr, typedict): - return self.infer_for_children(expr, typedict, - children=[expr.condition, expr.then, expr.else_]) - - def map_comparison(self, expr, typedict): - return self.infer_for_children(expr, typedict, - children=[expr.left, expr.right]) - - def map_if_positive(self, expr, typedict): - return self.infer_for_children(expr, typedict, - children=[expr.criterion, expr.then, expr.else_]) - - def map_call(self, expr, typedict): - # assumes functions to be non-type-changing - return self.infer_for_children(expr, typedict, - children=expr.parameters) - - def map_operator_binding(self, expr, typedict): - 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, - (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, - (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, - (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, ( - 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, sym.ElementwiseMaxOperator): - typedict[expr.field] = typedict[expr].unify( - type_info.KnownVolume(), expr.field) - return self.rec(expr.field, typedict) - - elif isinstance(expr.op, sym.RestrictToBoundary): - # upward propagation: argument has same rep tag as result - typedict[expr.field] = type_info.KnownVolume().unify( - extract_representation(type_info), expr.field) - - self.rec(expr.field, typedict) - - # downward propagation: result has same rep tag as argument - return type_info.KnownBoundary(expr.op.tag).unify( - extract_representation(typedict[expr.field]), expr) - - elif isinstance(expr.op, sym.FluxExchangeOperator): - raise NotImplementedError - - elif isinstance(expr.op, sym.FluxOperatorBase): - from pytools.obj_array import with_object_array_or_scalar - - repr_tag_cell = [type_info.no_type] - - def process_vol_flux_arg(flux_arg): - typedict[flux_arg] = type_info.KnownInteriorFaces() \ - .unify(repr_tag_cell[0], flux_arg) - repr_tag_cell[0] = extract_representation( - self.rec(flux_arg, typedict)) - - 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) - - repr_tag_cell[0] = extract_representation( - self.rec(flux_arg, typedict)) - - bpair = expr.field - with_object_array_or_scalar(process_vol_flux_arg, bpair.field) - with_object_array_or_scalar(process_bdry_flux_arg, bpair.bfield) - else: - with_object_array_or_scalar(process_vol_flux_arg, expr.field) - - return type_info.VolumeVector(NodalRepresentation()) - - 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, 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, sym.QuadratureBoundaryGridUpsampler): - typedict[expr.field] = type_info.BoundaryVector( - expr.op.boundary_tag, NodalRepresentation()) - self.rec(expr.field, typedict) - return type_info.BoundaryVector( - expr.op.boundary_tag, - QuadratureRepresentation(expr.op.quadrature_tag)) - - elif isinstance(expr.op, sym.ElementwiseLinearOperator): - typedict[expr.field] = type_info.VolumeVector(NodalRepresentation()) - self.rec(expr.field, typedict) - return type_info.VolumeVector(NodalRepresentation()) - - else: - raise RuntimeError("TypeInferrer doesn't know how to handle '%s'" - % expr.op) - - def map_whole_domain_flux(self, expr, typedict): - repr_tag_cell = [type_info.no_type] - - def process_vol_flux_arg(flux_arg): - typedict[flux_arg] = type_info.KnownInteriorFaces() \ - .unify(repr_tag_cell[0], flux_arg) - repr_tag_cell[0] = extract_representation( - self.rec(flux_arg, typedict)) - - def process_bdry_flux_arg(flux_arg): - typedict[flux_arg] = type_info.KnownBoundary(bpair.tag) \ - .unify(repr_tag_cell[0], flux_arg) - - repr_tag_cell[0] = extract_representation( - self.rec(flux_arg, typedict)) - - from pytools.obj_array import with_object_array_or_scalar - for int_flux_info in expr.interiors: - with_object_array_or_scalar(process_vol_flux_arg, - int_flux_info.field_expr) - - for bdry_flux_info in expr.boundaries: - bpair = bdry_flux_info.bpair - with_object_array_or_scalar(process_vol_flux_arg, bpair.field) - with_object_array_or_scalar(process_bdry_flux_arg, bpair.bfield) - - return type_info.VolumeVector(NodalRepresentation()) - - def map_flux_exchange(self, expr, typedict): - for arg in expr.arg_fields: - typedict[arg] = type_info.VolumeVector(NodalRepresentation()) - - from grudge.mesh import TAG_RANK_BOUNDARY - return type_info.BoundaryVector( - TAG_RANK_BOUNDARY(expr.rank), - NodalRepresentation()) - - def map_constant(self, expr, typedict): - return type_info.Scalar().unify(typedict[expr], expr) - - def map_variable(self, expr, typedict): - # user-facing variables are nodal - return type_info.KnownRepresentation(NodalRepresentation())\ - .unify(typedict[expr], expr) - - map_subscript = map_variable - - def map_scalar_parameter(self, expr, typedict): - return type_info.Scalar().unify(typedict[expr], expr) - - def map_ones(self, expr, typedict): - # FIXME: This is a bit dumb. If the quadrature_tag is None, - # we don't know whether the expression was specialized - # to 'no quadrature' or if it simply does not know yet - # whether it will be on a quadrature grid. - if expr.quadrature_tag is not None: - return (type_info.VolumeVector( - QuadratureRepresentation(expr.quadrature_tag)) - .unify(typedict[expr], expr)) - else: - return (type_info.VolumeVector(NodalRepresentation()) - .unify(typedict[expr], expr)) - - def map_node_coordinate_component(self, expr, typedict): - # FIXME: This is a bit dumb. If the quadrature_tag is None, - # we don't know whether the expression was specialized - # to 'no quadrature' or if it simply does not know yet - # whether it will be on a quadrature grid. - if expr.quadrature_tag is not None: - return (type_info.VolumeVector( - QuadratureRepresentation(expr.quadrature_tag)) - .unify(typedict[expr], expr)) - else: - return (type_info.VolumeVector(NodalRepresentation()) - .unify(typedict[expr], expr)) - - def map_normal_component(self, expr, typedict): - # FIXME: This is a bit dumb. If the quadrature_tag is None, - # we don't know whether the expression was specialized - # to 'no quadrature' or if it simply does not know yet - # whether it will be on a quadrature grid. - - if expr.quadrature_tag is not None: - return (type_info.BoundaryVector(expr.boundary_tag, - QuadratureRepresentation(expr.quadrature_tag)) - .unify(typedict[expr], expr)) - else: - return (type_info.KnownBoundary(expr.boundary_tag) - .unify(typedict[expr], expr)) - - def map_jacobian(self, expr, typedict): - return type_info.KnownVolume() - - map_forward_metric_derivative = map_jacobian - map_inverse_metric_derivative = map_jacobian - - def map_common_subexpression(self, expr, typedict): - outer_tp = typedict[expr] - - last_tp = self.cse_last_results.get(expr, type_info.no_type) - if outer_tp != last_tp or last_tp == type_info.no_type: - # re-run inner type inference with new outer information - typedict[expr.child] = outer_tp - new_tp = self.rec(expr.child, typedict) - - # For correct caching, we need to make sure that - # information below this level has fully propagated. - while True: - typedict.change_flag = False - new_tp = self.rec(expr.child, typedict) - - if not typedict.change_flag: - # nothing has changed any more - break - - self.cse_last_results[expr] = new_tp - typedict[expr.child] = new_tp - - # we can be sure we *have* changed something - typedict.change_flag = True - return new_tp - else: - return last_tp - -# }}} - - -# vim: foldmethod=marker diff --git a/grudge/symbolic/primitives.py b/grudge/symbolic/primitives.py index 0db0a8908a2b9647c87e15046480de8d3904fbe5..11d5ae8a5dc802b0a75d0f705494b6483ac40b78 100644 --- a/grudge/symbolic/primitives.py +++ b/grudge/symbolic/primitives.py @@ -73,9 +73,9 @@ Symbols ^^^^^^^ .. autoclass:: Variable +.. autoclass:: ScalarVariable .. autoclass:: make_sym_array .. autoclass:: make_sym_mv -.. autoclass:: ScalarParameter .. autoclass:: CFunction Helpers @@ -191,6 +191,9 @@ class DOFDesc(object): if domain_tag is DTAG_SCALAR and quadrature_tag is not None: raise ValueError("cannot have nontrivial quadrature tag on scalar") + if quadrature_tag is None: + quadrature_tag = QTAG_NONE + self.domain_tag = domain_tag self.quadrature_tag = quadrature_tag @@ -297,6 +300,9 @@ class cse_scope(cse_scope_base): # noqa class Variable(HasDOFDesc, ExpressionBase, pymbolic.primitives.Variable): + """A user-supplied input variable with a known :class:`DOFDesc`. + """ + def __init__(self, name, dd=None): if dd is None: dd = DD_VOLUME diff --git a/test/test_grudge.py b/test/test_grudge.py index 3372567408a29c3a3927619df8f6f5bf9ecd14b7..0ceaf8fec8339734e93dbdbfd12bdd1e603f171f 100644 --- a/test/test_grudge.py +++ b/test/test_grudge.py @@ -106,9 +106,9 @@ def test_1d_mass_mat_trig(ctx_factory): mass_op = bind(discr, sym.MassOperator()(sym.var("f"))) - num_integral_1 = np.dot(ones.get(), mass_op(queue, f=f).get()) - num_integral_2 = np.dot(f.get(), mass_op(queue, f=ones).get()) - num_integral_3 = bind(discr, sym.integral(sym.var("f")))(queue, f=f).get() + num_integral_1 = np.dot(ones.get(), mass_op(queue, f=f)) + num_integral_2 = np.dot(f.get(), mass_op(queue, f=ones)) + num_integral_3 = bind(discr, sym.integral(sym.var("f")))(queue, f=f) true_integral = 13*np.pi/2 err_1 = abs(num_integral_1-true_integral) @@ -200,7 +200,7 @@ def test_2d_gauss_theorem(ctx_factory): dd=sym.BTAG_ALL) )(queue) - assert abs(gauss_err.get()) < 1e-13 + assert abs(gauss_err) < 1e-13 @pytest.mark.parametrize(("mesh_name", "mesh_pars"), [ @@ -211,6 +211,7 @@ def test_2d_gauss_theorem(ctx_factory): @pytest.mark.parametrize("op_type", ["strong", "weak"]) @pytest.mark.parametrize("flux_type", ["upwind"]) @pytest.mark.parametrize("order", [3, 4, 5]) +# test: 'test_convergence_advec(cl._csc, "disk", [0.1, 0.05], "strong", "upwind", 3)' def test_convergence_advec(ctx_factory, mesh_name, mesh_pars, op_type, flux_type, order, visualize=False): """Test whether 2D advection actually converges""" @@ -322,7 +323,7 @@ def test_convergence_advec(ctx_factory, mesh_name, mesh_pars, op_type, flux_type error_l2 = bind(discr, sym.norm(2, sym.var("u")-u_analytic(sym.nodes(dim))))( - queue, t=last_t, u=last_u).get() + queue, t=last_t, u=last_u) print(h, error_l2) eoc_rec.add_data_point(h, error_l2)