From 706e6753f90fcc42c712f17b94e7306b5073d2d6 Mon Sep 17 00:00:00 2001
From: Andreas Kloeckner <inform@tiker.net>
Date: Fri, 10 Jun 2022 12:52:38 -0500
Subject: [PATCH] Drop symbolic execution in grudge
---
examples/old_symbolics/dagrt-fusion.py | 1447 --------------------
examples/old_symbolics/symbolic_wave_op.py | 166 ---
grudge/__init__.py | 5 +-
grudge/discretization.py | 11 -
grudge/execution.py | 822 -----------
grudge/function_registry.py | 193 ---
grudge/symbolic/__init__.py | 29 -
grudge/symbolic/compiler.py | 1369 ------------------
grudge/symbolic/dofdesc_inference.py | 232 ----
grudge/symbolic/mappers/__init__.py | 1298 ------------------
grudge/symbolic/operators.py | 862 ------------
grudge/symbolic/primitives.py | 724 ----------
grudge/symbolic/tools.py | 141 --
grudge/trace_pair.py | 5 -
test/test_grudge.py | 3 +-
test/test_grudge_sym_old.py | 943 -------------
16 files changed, 3 insertions(+), 8247 deletions(-)
delete mode 100755 examples/old_symbolics/dagrt-fusion.py
delete mode 100644 examples/old_symbolics/symbolic_wave_op.py
delete mode 100644 grudge/execution.py
delete mode 100644 grudge/function_registry.py
delete mode 100644 grudge/symbolic/__init__.py
delete mode 100644 grudge/symbolic/compiler.py
delete mode 100644 grudge/symbolic/dofdesc_inference.py
delete mode 100644 grudge/symbolic/mappers/__init__.py
delete mode 100644 grudge/symbolic/operators.py
delete mode 100644 grudge/symbolic/primitives.py
delete mode 100644 grudge/symbolic/tools.py
delete mode 100644 test/test_grudge_sym_old.py
diff --git a/examples/old_symbolics/dagrt-fusion.py b/examples/old_symbolics/dagrt-fusion.py
deleted file mode 100755
index 6ba48031..00000000
--- a/examples/old_symbolics/dagrt-fusion.py
+++ /dev/null
@@ -1,1447 +0,0 @@
-#!/usr/bin/env python3
-"""Study of operator fusion (inlining) for time integration operators in Grudge.
-"""
-
-__copyright__ = """
-Copyright (C) 2015 Andreas Kloeckner
-Copyright (C) 2019 Matt Wala
-"""
-
-__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.
-"""
-
-# FIXME:
-# Results before https://github.com/inducer/grudge/pull/15 were better:
-#
-# Operator | \parbox{1in}{\centering \% Memory Ops. Due to Scalar Assignments}
-# -------------+-------------------------------------------------------------------
-# 2D: Baseline | 51.1
-# 2D: Inlined | 48.9
-# 3D: Baseline | 50.1
-# 3D: Inlined | 48.6
-# INFO:__main__:Wrote '<stdout>'
-# ==== Scalar Assigment Inlining Impact ====
-# Operator | Bytes Read | Bytes Written | Total | \% of Baseline
-# -------------+------------+---------------+------------+----------------
-# 2D: Baseline | 9489600 | 3348000 | 12837600 | 100
-# 2D: Inlined | 8949600 | 2808000 | 11757600 | 91.6
-# 3D: Baseline | 1745280000 | 505440000 | 2250720000 | 100
-# 3D: Inlined | 1680480000 | 440640000 | 2121120000 | 94.2
-# INFO:__main__:Wrote '<stdout>'
-
-
-import contextlib
-import logging
-import numpy as np
-import os
-import sys
-import pyopencl as cl
-import pyopencl.array # noqa
-import pyopencl.tools as cl_tools
-import pytest
-
-import dagrt.language as lang
-import pymbolic.primitives as p
-from pymbolic.mapper import IdentityMapper
-
-from arraycontext import PyOpenCLArrayContext
-
-from meshmode.dof_array import DOFArray
-
-import grudge.dof_desc as dof_desc
-import grudge.symbolic.mappers as gmap
-import grudge.symbolic.operators as sym_op
-
-from grudge.execution import ExecutionMapper
-from grudge.function_registry import base_function_registry
-from pymbolic.mapper import Mapper
-from pymbolic.mapper.evaluator import EvaluationMapper \
- as PymbolicEvaluationMapper
-from pytools import memoize
-from pytools.obj_array import flat_obj_array
-
-from grudge import sym, bind, DiscretizationCollection
-from leap.rk import LSRK4MethodBuilder
-
-from pyopencl.tools import ( # noqa
- pytest_generate_tests_for_pyopencl as pytest_generate_tests)
-
-
-logging.basicConfig(level=logging.INFO)
-
-logger = logging.getLogger(__name__)
-
-
-SKIP_TESTS = int(os.environ.get("SKIP_TESTS", 0))
-PAPER_OUTPUT = int(os.environ.get("PAPER_OUTPUT", 0))
-OUT_DIR = os.environ.get("OUT_DIR", ".")
-
-
-@contextlib.contextmanager
-def open_output_file(filename):
- if not PAPER_OUTPUT:
- yield sys.stdout
- sys.stdout.flush()
- else:
- try:
- outfile = open(os.path.join(OUT_DIR, filename), "w")
- yield outfile
- finally:
- outfile.close()
-
-
-def dof_array_nbytes(ary: np.ndarray):
- if isinstance(ary, np.ndarray) and ary.dtype.char == "O":
- return sum(
- dof_array_nbytes(ary[idx])
- for idx in np.ndindex(ary.shape))
- elif isinstance(ary, DOFArray):
- return sum(dof_array_nbytes(ary_i) for ary_i in ary)
- else:
- return ary.nbytes
-
-
-# {{{ topological sort
-
-def topological_sort(stmts, root_deps):
- id_to_stmt = {stmt.id: stmt for stmt in stmts}
-
- ordered_stmts = []
- satisfied = set()
-
- def satisfy_dep(name):
- if name in satisfied:
- return
-
- stmt = id_to_stmt[name]
- for dep in sorted(stmt.depends_on):
- satisfy_dep(dep)
- ordered_stmts.append(stmt)
- satisfied.add(name)
-
- for d in root_deps:
- satisfy_dep(d)
-
- return ordered_stmts
-
-# }}}
-
-
-# {{{ leap to grudge translation
-
-# Use evaluation, not identity mappers to propagate symbolic vectors to
-# outermost level.
-
-class DagrtToGrudgeRewriter(PymbolicEvaluationMapper):
- def __init__(self, context):
- self.context = context
-
- def map_variable(self, expr):
- return self.context[expr.name]
-
- def map_call(self, expr):
- raise ValueError("function call not expected")
-
-
-class GrudgeArgSubstitutor(gmap.SymbolicEvaluator):
- def __init__(self, args):
- super().__init__(context={})
- self.args = args
-
- def map_grudge_variable(self, expr):
- if expr.name in self.args:
- return self.args[expr.name]
- return super().map_variable(expr)
-
-
-# {{{ isolate function calls
-
-# (copied from dagrt pre-https://github.com/inducer/dagrt/pull/12)
-
-class FunctionCallIsolator(IdentityMapper):
- def __init__(self, new_statements,
- stmt_id_gen, var_name_gen):
- super().__init__()
- self.new_statements = new_statements
- self.stmt_id_gen = stmt_id_gen
- self.var_name_gen = var_name_gen
-
- def isolate_call(self, expr, base_condition, base_deps, extra_deps,
- super_method):
- # FIXME: These aren't awesome identifiers.
- tmp_var_name = self.var_name_gen("tmp")
-
- tmp_stmt_id = self.stmt_id_gen("tmp")
- extra_deps.append(tmp_stmt_id)
-
- sub_extra_deps = []
- rec_result = super_method(
- expr, base_deps, sub_extra_deps)
-
- from pymbolic.primitives import Call, CallWithKwargs
- assert isinstance(rec_result, (Call, CallWithKwargs))
-
- parameters = []
- kw_parameters = {}
-
- for par in rec_result.parameters:
- parameters.append(par)
-
- if isinstance(rec_result, CallWithKwargs):
- for par_name, par in rec_result.kw_parameters.items():
- kw_parameters[par_name] = par
-
- from dagrt.language import AssignFunctionCall
- new_stmt = AssignFunctionCall(
- assignees=(tmp_var_name,),
- function_id=rec_result.function.name,
- parameters=tuple(parameters),
- kw_parameters=kw_parameters,
- id=tmp_stmt_id,
- condition=base_condition,
- depends_on=base_deps | frozenset(sub_extra_deps))
-
- self.new_statements.append(new_stmt)
-
- from pymbolic import var
- return var(tmp_var_name)
-
- def map_call(self, expr, base_condition, base_deps, extra_deps):
- return self.isolate_call(
- expr, base_condition, base_deps, extra_deps,
- super().map_call)
-
- def map_call_with_kwargs(self, expr, base_condition, base_deps, extra_deps):
- return self.isolate_call(
- expr, base_condition, base_deps, extra_deps,
- super()
- .map_call_with_kwargs)
-
-
-def isolate_function_calls_in_phase(phase, stmt_id_gen, var_name_gen):
- new_statements = []
-
- fci = FunctionCallIsolator(
- new_statements=new_statements,
- stmt_id_gen=stmt_id_gen,
- var_name_gen=var_name_gen)
-
- for stmt in sorted(phase.statements, key=lambda stmt_: stmt_.id):
- new_deps = []
-
- from dagrt.language import Assign
- if isinstance(stmt, Assign):
- assert not stmt.loops
- new_statements.append(
- stmt
- .map_expressions(
- lambda expr: fci(
- expr, stmt.condition, stmt.depends_on, new_deps))
- .copy(depends_on=stmt.depends_on | frozenset(new_deps)))
- from pymbolic.primitives import Call, CallWithKwargs
- assert not isinstance(new_statements[-1].rhs,
- (Call, CallWithKwargs))
- else:
- new_statements.append(stmt)
-
- return phase.copy(statements=new_statements)
-
-
-def isolate_function_calls(dag):
- """
- :func:`isolate_function_arguments` should be
- called before this.
- """
-
- stmt_id_gen = dag.get_stmt_id_generator()
- var_name_gen = dag.get_var_name_generator()
-
- new_phases = {}
- for phase_name, phase in dag.phases.items():
- new_phases[phase_name] = isolate_function_calls_in_phase(
- phase, stmt_id_gen, var_name_gen)
-
- return dag.copy(phases=new_phases)
-
-# }}}
-
-
-def transcribe_phase(dag, field_var_name, field_components, phase_name,
- sym_operator):
- """Generate a Grudge operator for a Dagrt time integrator phase.
-
- Arguments:
-
- dag: The Dagrt code object for the time integrator
-
- field_var_name: The name of the simulation variable
-
- field_components: The number of components (fields) in the variable
-
- phase_name: The name of the phase to transcribe
-
- sym_operator: The Grudge symbolic expression to substitue for the
- right-hand side evaluation in the Dagrt code
- """
- sym_operator = gmap.OperatorBinder()(sym_operator)
- phase = dag.phases[phase_name]
-
- ctx = {
- "<t>": sym.var("input_t", dof_desc.DD_SCALAR),
- "<dt>": sym.var("input_dt", dof_desc.DD_SCALAR),
- f"<state>{field_var_name}": sym.make_sym_array(
- f"input_{field_var_name}", field_components),
- "<p>residual": sym.make_sym_array(
- "input_residual", field_components),
- }
-
- rhs_name = f"<func>{field_var_name}"
- output_vars = [v for v in ctx]
- yielded_states = []
-
- ordered_stmts = topological_sort(
- isolate_function_calls_in_phase(
- phase,
- dag.get_stmt_id_generator(),
- dag.get_var_name_generator()).statements,
- phase.depends_on)
-
- for stmt in ordered_stmts:
- if stmt.condition is not True:
- raise NotImplementedError(
- "non-True condition (in statement '%s') not supported"
- % stmt.id)
-
- if isinstance(stmt, lang.Nop):
- pass
-
- elif isinstance(stmt, lang.Assign):
- if not isinstance(stmt.lhs, p.Variable):
- raise NotImplementedError("lhs of statement %s is not a variable: %s"
- % (stmt.id, stmt.lhs))
- ctx[stmt.lhs.name] = sym.cse(
- DagrtToGrudgeRewriter(ctx)(stmt.rhs),
- (
- stmt.lhs.name
- .replace("<", "")
- .replace(">", "")))
-
- elif isinstance(stmt, lang.AssignFunctionCall):
- if stmt.function_id != rhs_name:
- raise NotImplementedError(
- "statement '%s' calls unsupported function '%s'"
- % (stmt.id, stmt.function_id))
-
- if stmt.parameters:
- raise NotImplementedError(
- "statement '%s' calls function '%s' with positional arguments"
- % (stmt.id, stmt.function_id))
-
- kwargs = {name: sym.cse(DagrtToGrudgeRewriter(ctx)(arg))
- for name, arg in stmt.kw_parameters.items()}
-
- if len(stmt.assignees) != 1:
- raise NotImplementedError(
- "statement '%s' calls function '%s' "
- "with more than one LHS"
- % (stmt.id, stmt.function_id))
-
- assignee, = stmt.assignees
- ctx[assignee] = GrudgeArgSubstitutor(kwargs)(sym_operator)
-
- elif isinstance(stmt, lang.YieldState):
- d2g = DagrtToGrudgeRewriter(ctx)
- yielded_states.append(
- (
- stmt.time_id,
- d2g(stmt.time),
- stmt.component_id,
- d2g(stmt.expression)))
-
- else:
- raise NotImplementedError("statement %s is of unsupported type ''%s'"
- % (stmt.id, type(stmt).__name__))
-
- return output_vars, [ctx[ov] for ov in output_vars], yielded_states
-
-# }}}
-
-
-# {{{ time integrator implementations
-
-class RK4TimeStepperBase:
-
- def __init__(self, component_getter):
- self.component_getter = component_getter
-
- def get_initial_context(self, fields, t_start, dt):
-
- # Flatten fields.
- flattened_fields = []
- for field in fields:
- if isinstance(field, list):
- flattened_fields.extend(field)
- else:
- flattened_fields.append(field)
- flattened_fields = flat_obj_array(*flattened_fields)
- del fields
-
- return {
- "input_t": t_start,
- "input_dt": dt,
- self.state_name: flattened_fields,
- "input_residual": flattened_fields,
- }
-
- def set_up_stepper(self, discr, field_var_name, sym_rhs, num_fields,
- function_registry=base_function_registry,
- exec_mapper_factory=ExecutionMapper):
- dt_method = LSRK4MethodBuilder(component_id=field_var_name)
- dt_code = dt_method.generate()
- self.field_var_name = field_var_name
- self.state_name = f"input_{field_var_name}"
-
- # Transcribe the phase.
- output_vars, results, yielded_states = transcribe_phase(
- dt_code, field_var_name, num_fields,
- "primary", sym_rhs)
-
- # Build the bound operator for the time integrator.
- output_t = results[0]
- output_dt = results[1]
- output_states = results[2]
- output_residuals = results[3]
-
- assert len(output_states) == num_fields
- assert len(output_states) == len(output_residuals)
-
- flattened_results = flat_obj_array(output_t, output_dt, *output_states)
-
- self.bound_op = bind(
- discr, flattened_results,
- function_registry=function_registry,
- exec_mapper_factory=exec_mapper_factory)
-
- def run(self, fields, t_start, dt, t_end, return_profile_data=False):
- context = self.get_initial_context(fields, t_start, dt)
-
- t = t_start
-
- while t <= t_end:
- if return_profile_data:
- profile_data = dict()
- else:
- profile_data = None
-
- results = self.bound_op(
- profile_data=profile_data,
- **context)
-
- if return_profile_data:
- results = results[0]
-
- t = results[0]
- context["input_t"] = t
- context["input_dt"] = results[1]
- output_states = results[2:]
- context[self.state_name] = output_states
-
- result = (t, self.component_getter(output_states))
- if return_profile_data:
- result += (profile_data,)
-
- yield result
-
-
-class RK4TimeStepper(RK4TimeStepperBase):
-
- def __init__(self, discr, field_var_name, grudge_bound_op,
- num_fields, component_getter, exec_mapper_factory=ExecutionMapper):
- """Arguments:
-
- field_var_name: The name of the simulation variable
-
- grudge_bound_op: The BoundExpression for the right-hand side
-
- num_fields: The number of components in the simulation variable
-
- component_getter: A function, which, given an object array
- representing the simulation variable, splits the array into
- its components
-
- """
- super().__init__(component_getter)
-
- # Construct sym_rhs to have the effect of replacing the RHS calls in the
- # dagrt code with calls of the grudge operator.
- from grudge.symbolic.primitives import FunctionSymbol, Variable
- call = sym.cse(
- FunctionSymbol("grudge_op")(*(
- (Variable("t", dd=dof_desc.DD_SCALAR),)
- + tuple(
- Variable(field_var_name, dd=dof_desc.DD_VOLUME)[i]
- for i in range(num_fields)))))
- sym_rhs = flat_obj_array(*(call[i] for i in range(num_fields)))
-
- self.grudge_bound_op = grudge_bound_op
-
- from grudge.function_registry import register_external_function
-
- freg = register_external_function(
- base_function_registry,
- "grudge_op",
- implementation=self._bound_op,
- dd=dof_desc.DD_VOLUME)
-
- self.set_up_stepper(
- discr, field_var_name, sym_rhs, num_fields,
- freg,
- exec_mapper_factory)
-
- self.component_getter = component_getter
-
- def _bound_op(self, array_context, t, *args, profile_data=None):
- context = {
- "t": t,
- self.field_var_name: flat_obj_array(*args)}
- result = self.grudge_bound_op(
- array_context, profile_data=profile_data, **context)
- if profile_data is not None:
- result = result[0]
- return result
-
- def get_initial_context(self, fields, t_start, dt):
- context = super().get_initial_context(fields, t_start, dt)
- context["grudge_op"] = self._bound_op
- return context
-
-
-class FusedRK4TimeStepper(RK4TimeStepperBase):
-
- def __init__(self, discr, field_var_name, sym_rhs, num_fields,
- component_getter, exec_mapper_factory=ExecutionMapper):
- super().__init__(component_getter)
- self.set_up_stepper(
- discr, field_var_name, sym_rhs, num_fields,
- base_function_registry,
- exec_mapper_factory)
-
-# }}}
-
-
-# {{{ problem setup code
-
-def _get_source_term(dims):
- source_center = np.array([0.1, 0.22, 0.33])[:dims]
- source_width = 0.05
- source_omega = 3
-
- sym_x = sym.nodes(dims)
- sym_source_center_dist = sym_x - source_center
- sym_t = sym.ScalarVariable("t")
-
- return (
- sym.sin(source_omega*sym_t)
- * sym.exp(
- -np.dot(sym_source_center_dist, sym_source_center_dist)
- / source_width**2))
-
-
-def get_wave_op_with_discr(actx, dims=2, order=4):
- from meshmode.mesh.generation import generate_regular_rect_mesh
- mesh = generate_regular_rect_mesh(
- a=(-0.5,)*dims,
- b=(0.5,)*dims,
- nelements_per_axis=(16,)*dims)
-
- logger.debug("%d elements", mesh.nelements)
-
- discr = DiscretizationCollection(actx, mesh, order=order)
-
- from symbolic_wave_op import WeakWaveOperator
- from meshmode.mesh import BTAG_ALL, BTAG_NONE
- op = WeakWaveOperator(0.1, dims,
- source_f=_get_source_term(dims),
- dirichlet_tag=BTAG_NONE,
- neumann_tag=BTAG_NONE,
- radiation_tag=BTAG_ALL,
- flux_type="upwind")
-
- op.check_bc_coverage(mesh)
-
- return (op.sym_operator(), discr)
-
-
-def get_wave_component(state_component):
- return (state_component[0], state_component[1:])
-
-# }}}
-
-
-# {{{ equivalence check between fused and non-fused versions
-
-def test_stepper_equivalence(ctx_factory, order=4):
- cl_ctx = ctx_factory()
- queue = cl.CommandQueue(cl_ctx)
- actx = PyOpenCLArrayContext(
- queue,
- allocator=cl_tools.MemoryPool(cl_tools.ImmediateAllocator(queue))
- )
-
- dims = 2
-
- sym_operator, discr = get_wave_op_with_discr(
- actx, dims=dims, order=order)
- #sym_operator_direct, discr = get_wave_op_with_discr_direct(
- # actx, dims=dims, order=order)
-
- if dims == 2:
- dt = 0.04
- elif dims == 3:
- dt = 0.02
-
- ic = flat_obj_array(discr.zeros(actx),
- [discr.zeros(actx) for i in range(discr.dim)])
-
- bound_op = bind(discr, sym_operator)
-
- stepper = RK4TimeStepper(
- discr, "w", bound_op, 1 + discr.dim, get_wave_component)
-
- fused_stepper = FusedRK4TimeStepper(
- discr, "w", sym_operator, 1 + discr.dim,
- get_wave_component)
-
- t_start = 0
- t_end = 0.5
- nsteps = int(np.ceil((t_end + 1e-9) / dt))
- print("dt=%g nsteps=%d" % (dt, nsteps))
-
- step = 0
-
- norm = bind(discr, sym.norm(2, sym.var("u_ref") - sym.var("u")))
-
- fused_steps = fused_stepper.run(ic, t_start, dt, t_end)
-
- for t_ref, (u_ref, _v_ref) in stepper.run(ic, t_start, dt, t_end):
- step += 1
- logger.debug("step %d/%d", step, nsteps)
- t, (u, v) = next(fused_steps)
- assert t == t_ref, step
- assert norm(u=u, u_ref=u_ref) <= 1e-13, step
-
-# }}}
-
-
-# {{{ execution mapper wrapper
-
-class ExecutionMapperWrapper(Mapper):
-
- def __init__(self, array_context, context, bound_op):
- self.inner_mapper = ExecutionMapper(array_context, context, bound_op)
- self.array_context = array_context
- self.context = context
- self.bound_op = bound_op
-
- def map_variable(self, expr):
- # Needed, because bound op execution can ask for variable values.
- return self.inner_mapper.map_variable(expr)
-
- def map_node_coordinate_component(self, expr):
- return self.inner_mapper.map_node_coordinate_component(expr)
-
- def map_grudge_variable(self, expr):
- # See map_variable()
- return self.inner_mapper.map_grudge_variable(expr)
-
-# }}}
-
-
-# {{{ mem op counter implementation
-
-class ExecutionMapperWithMemOpCounting(ExecutionMapperWrapper):
- # This is a skeleton implementation that only has just enough functionality
- # for the wave-min example to work.
-
- # {{{ expressions
-
- def map_profiled_call(self, expr, profile_data):
- args = [self.inner_mapper.rec(p) for p in expr.parameters]
- return self.inner_mapper.function_registry[expr.function.name](
- self.array_context, *args, profile_data=profile_data)
-
- def map_profiled_essentially_elementwise_linear(self, op, field_expr,
- profile_data):
- result = getattr(self.inner_mapper, op.mapper_method)(op, field_expr)
-
- if profile_data is not None:
- # We model the cost to load the input and write the output. In
- # particular, we assume the elementwise matrices are negligible in
- # size and thus ignorable.
-
- field = self.inner_mapper.rec(field_expr)
- profile_data["bytes_read"] = (
- profile_data.get("bytes_read", 0)
- + dof_array_nbytes(field))
- profile_data["bytes_written"] = (
- profile_data.get("bytes_written", 0)
- + dof_array_nbytes(result))
-
- if op.mapper_method == "map_projection":
- profile_data["interp_bytes_read"] = (
- profile_data.get("interp_bytes_read", 0)
- + dof_array_nbytes(field))
- profile_data["interp_bytes_written"] = (
- profile_data.get("interp_bytes_written", 0)
- + dof_array_nbytes(result))
-
- return result
-
- # }}}
-
- # {{{ instruction mappings
-
- def process_assignment_expr(self, expr, profile_data):
- if isinstance(expr, p.Call):
- assert expr.mapper_method == "map_call"
- val = self.map_profiled_call(expr, profile_data)
-
- elif isinstance(expr, sym.OperatorBinding):
- if isinstance(
- expr.op,
- (
- # TODO: Not comprehensive.
- sym_op.ProjectionOperator,
- sym_op.RefFaceMassOperator,
- sym_op.RefMassOperator,
- sym_op.RefInverseMassOperator,
- sym_op.OppositeInteriorFaceSwap)):
- val = self.map_profiled_essentially_elementwise_linear(
- expr.op, expr.field, profile_data)
-
- else:
- raise AssertionError("unknown operator: %s" % expr.op)
-
- else:
- logger.debug("assignment not profiled: %s", expr)
- val = self.inner_mapper.rec(expr)
-
- return val
-
- def map_insn_assign(self, insn, profile_data):
- result = []
- for name, expr in zip(insn.names, insn.exprs):
- result.append((name, self.process_assignment_expr(expr, profile_data)))
- return result, []
-
- def map_insn_loopy_kernel(self, insn, profile_data):
- kdescr = insn.kernel_descriptor
- discr = self.inner_mapper.dcoll.discr_from_dd(kdescr.governing_dd)
-
- dof_array_kwargs = {}
- other_kwargs = {}
-
- for name, expr in kdescr.input_mappings.items():
- v = self.inner_mapper.rec(expr)
- if isinstance(v, DOFArray):
- dof_array_kwargs[name] = v
-
- if profile_data is not None:
- size = dof_array_nbytes(v)
- profile_data["bytes_read"] = (
- profile_data.get("bytes_read", 0) + size)
- profile_data["bytes_read_by_scalar_assignments"] = (
- profile_data.get("bytes_read_by_scalar_assignments", 0)
- + size)
- else:
- assert not isinstance(v, np.ndarray)
- other_kwargs[name] = v
-
- for name in kdescr.scalar_args():
- v = other_kwargs[name]
- if isinstance(v, (int, float)):
- other_kwargs[name] = discr.real_dtype.type(v)
- elif isinstance(v, complex):
- other_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)))
-
- result = {}
- for grp in discr.groups:
- kwargs = other_kwargs.copy()
- kwargs["nelements"] = grp.nelements
- kwargs["nunit_dofs"] = grp.nunit_dofs
-
- for name, ary in dof_array_kwargs.items():
- kwargs[name] = ary[grp.index]
-
- knl_result = self.inner_mapper.array_context.call_loopy(
- kdescr.loopy_kernel, **kwargs)
-
- for name, val in knl_result.items():
- result.setdefault(name, []).append(val)
-
- result = {
- name: DOFArray(self.inner_mapper.array_context, tuple(val))
- for name, val in result.items()}
-
- for val in result.values():
- assert isinstance(val, DOFArray)
- if profile_data is not None:
- size = dof_array_nbytes(val)
- profile_data["bytes_written"] = (
- profile_data.get("bytes_written", 0) + size)
- profile_data["bytes_written_by_scalar_assignments"] = (
- profile_data.get("bytes_written_by_scalar_assignments", 0)
- + size)
-
- return list(result.items()), []
-
- def map_insn_assign_to_discr_scoped(self, insn, profile_data=None):
- assignments = []
-
- for name, expr in zip(insn.names, insn.exprs):
- logger.debug("assignment not profiled: %s <- %s", name, expr)
- inner_mapper = self.inner_mapper
- value = inner_mapper.rec(expr)
- inner_mapper.dcoll._discr_scoped_subexpr_name_to_value[name] = value
- assignments.append((name, value))
-
- return assignments, []
-
- def map_insn_assign_from_discr_scoped(self, insn, profile_data=None):
- return [(
- insn.name,
- self.inner_mapper.
- dcoll._discr_scoped_subexpr_name_to_value[insn.name])], []
-
- def map_insn_rank_data_swap(self, insn, profile_data):
- raise NotImplementedError("no profiling for instruction: %s" % insn)
-
- def map_insn_diff_batch_assign(self, insn, profile_data):
- assignments, futures = self.inner_mapper.map_insn_diff_batch_assign(insn)
-
- if profile_data is not None:
- # We model the cost to load the input and write the output. In
- # particular, we assume the elementwise matrices are negligible in
- # size and thus ignorable.
-
- field = self.inner_mapper.rec(insn.field)
- profile_data["bytes_read"] = (
- profile_data.get("bytes_read", 0) + dof_array_nbytes(field))
-
- for _, value in assignments:
- profile_data["bytes_written"] = (
- profile_data.get("bytes_written", 0)
- + dof_array_nbytes(value))
-
- return assignments, futures
-
- # }}}
-
-# }}}
-
-
-# {{{ mem op counter check
-
-def test_assignment_memory_model(ctx_factory):
- cl_ctx = ctx_factory()
- queue = cl.CommandQueue(cl_ctx)
- actx = PyOpenCLArrayContext(
- queue,
- allocator=cl_tools.MemoryPool(cl_tools.ImmediateAllocator(queue))
- )
-
- _, discr = get_wave_op_with_discr(actx, dims=2, order=3)
-
- # Assignment instruction
- bound_op = bind(
- discr,
- sym.Variable("input0", dof_desc.DD_VOLUME)
- + sym.Variable("input1", dof_desc.DD_VOLUME),
- exec_mapper_factory=ExecutionMapperWithMemOpCounting)
-
- input0 = discr.zeros(actx)
- input1 = discr.zeros(actx)
-
- result, profile_data = bound_op(
- profile_data={},
- input0=input0,
- input1=input1)
-
- assert profile_data["bytes_read"] == \
- dof_array_nbytes(input0) + dof_array_nbytes(input1)
- assert profile_data["bytes_written"] == dof_array_nbytes(result)
-
-
-@pytest.mark.parametrize("use_fusion", (True, False))
-def test_stepper_mem_ops(ctx_factory, use_fusion):
- cl_ctx = ctx_factory()
- queue = cl.CommandQueue(cl_ctx)
- actx = PyOpenCLArrayContext(
- queue,
- allocator=cl_tools.MemoryPool(cl_tools.ImmediateAllocator(queue))
- )
-
- dims = 2
-
- sym_operator, discr = get_wave_op_with_discr(
- actx, dims=dims, order=3)
-
- t_start = 0
- dt = 0.04
- t_end = 0.2
-
- ic = flat_obj_array(discr.zeros(actx),
- [discr.zeros(actx) for i in range(discr.dim)])
-
- if not use_fusion:
- bound_op = bind(
- discr, sym_operator,
- exec_mapper_factory=ExecutionMapperWithMemOpCounting)
-
- stepper = RK4TimeStepper(
- discr, "w", bound_op, 1 + discr.dim,
- get_wave_component,
- exec_mapper_factory=ExecutionMapperWithMemOpCounting)
-
- else:
- stepper = FusedRK4TimeStepper(
- discr, "w", sym_operator, 1 + discr.dim,
- get_wave_component,
- exec_mapper_factory=ExecutionMapperWithMemOpCounting)
-
- step = 0
-
- nsteps = int(np.ceil((t_end + 1e-9) / dt))
- for (_, _, profile_data) in stepper.run( # noqa: B007
- ic, t_start, dt, t_end, return_profile_data=True):
- step += 1
- logger.info("step %d/%d", step, nsteps)
-
- logger.info("using fusion? %s", use_fusion)
- logger.info("bytes read: %d", profile_data["bytes_read"])
- logger.info("bytes written: %d", profile_data["bytes_written"])
- logger.info("bytes total: %d",
- profile_data["bytes_read"] + profile_data["bytes_written"])
-
-# }}}
-
-
-# {{{ execution mapper with timing
-
-SECONDS_PER_NANOSECOND = 10**9
-
-
-class TimingFuture:
-
- def __init__(self, start_event, stop_event):
- self.start_event = start_event
- self.stop_event = stop_event
-
- def elapsed(self):
- cl.wait_for_events([self.start_event, self.stop_event])
- return (
- self.stop_event.profile.end
- - self.start_event.profile.end) / SECONDS_PER_NANOSECOND
-
-
-from collections.abc import MutableSequence
-
-
-class TimingFutureList(MutableSequence):
-
- def __init__(self, *args, **kwargs):
- self._list = list(*args, **kwargs)
-
- def __len__(self):
- return len(self._list)
-
- def __getitem__(self, idx):
- return self._list[idx]
-
- def __setitem__(self, idx, val):
- self._list[idx] = val
-
- def __delitem__(self, idx):
- del self._list[idx]
-
- def insert(self, idx, val):
- self._list.insert(idx, val)
-
- def elapsed(self):
- return sum(future.elapsed() for future in self._list)
-
-
-def time_insn(f):
- time_field_name = "time_%s" % f.__name__
-
- def wrapper(self, insn, profile_data):
- if profile_data is None:
- return f(self, insn, profile_data)
-
- start = cl.enqueue_marker(self.array_context.queue)
- retval = f(self, insn, profile_data)
- end = cl.enqueue_marker(self.array_context.queue)
- profile_data\
- .setdefault(time_field_name, TimingFutureList())\
- .append(TimingFuture(start, end))
-
- return retval
-
- return wrapper
-
-
-class ExecutionMapperWithTiming(ExecutionMapperWrapper):
-
- def map_profiled_call(self, expr, profile_data):
- args = [self.inner_mapper.rec(p) for p in expr.parameters]
- return self.inner_mapper.function_registry[expr.function.name](
- self.array_context, *args, profile_data=profile_data)
-
- def map_profiled_operator_binding(self, expr, profile_data):
- if profile_data is None:
- return self.inner_mapper.map_operator_binding(expr)
-
- start = cl.enqueue_marker(self.array_context.queue)
- retval = self.inner_mapper.map_operator_binding(expr)
- end = cl.enqueue_marker(self.array_context.queue)
- time_field_name = "time_op_%s" % expr.op.mapper_method
- profile_data\
- .setdefault(time_field_name, TimingFutureList())\
- .append(TimingFuture(start, end))
-
- return retval
-
- def map_insn_assign_to_discr_scoped(self, insn, profile_data):
- return self.inner_mapper.map_insn_assign_to_discr_scoped(insn, profile_data)
-
- def map_insn_assign_from_discr_scoped(self, insn, profile_data):
- return self.\
- inner_mapper.map_insn_assign_from_discr_scoped(insn, profile_data)
-
- @time_insn
- def map_insn_loopy_kernel(self, *args, **kwargs):
- return self.inner_mapper.map_insn_loopy_kernel(*args, **kwargs)
-
- def map_insn_assign(self, insn, profile_data):
- if len(insn.exprs) == 1:
- if isinstance(insn.exprs[0], p.Call):
- assert insn.exprs[0].mapper_method == "map_call"
- val = self.map_profiled_call(insn.exprs[0], profile_data)
- return [(insn.names[0], val)], []
- elif isinstance(insn.exprs[0], sym.OperatorBinding):
- assert insn.exprs[0].mapper_method == "map_operator_binding"
- val = self.map_profiled_operator_binding(insn.exprs[0], profile_data)
- return [(insn.names[0], val)], []
-
- return self.inner_mapper.map_insn_assign(insn, profile_data)
-
- @time_insn
- def map_insn_diff_batch_assign(self, insn, profile_data):
- return self.inner_mapper.map_insn_diff_batch_assign(insn, profile_data)
-
-# }}}
-
-
-# {{{ timing check
-
-@pytest.mark.parametrize("use_fusion", (True, False))
-def test_stepper_timing(ctx_factory, use_fusion):
- cl_ctx = ctx_factory()
- queue = cl.CommandQueue(
- cl_ctx,
- properties=cl.command_queue_properties.PROFILING_ENABLE)
- actx = PyOpenCLArrayContext(
- queue,
- allocator=cl_tools.MemoryPool(cl_tools.ImmediateAllocator(queue))
- )
-
- dims = 3
-
- sym_operator, discr = get_wave_op_with_discr(
- actx, dims=dims, order=3)
-
- t_start = 0
- dt = 0.04
- t_end = 0.1
-
- ic = flat_obj_array(discr.zeros(actx),
- [discr.zeros(actx) for i in range(discr.dim)])
-
- if not use_fusion:
- bound_op = bind(
- discr, sym_operator,
- exec_mapper_factory=ExecutionMapperWithTiming)
-
- stepper = RK4TimeStepper(
- discr, "w", bound_op, 1 + discr.dim,
- get_wave_component,
- exec_mapper_factory=ExecutionMapperWithTiming)
-
- else:
- stepper = FusedRK4TimeStepper(
- discr, "w", sym_operator, 1 + discr.dim,
- get_wave_component,
- exec_mapper_factory=ExecutionMapperWithTiming)
-
- step = 0
-
- import time
- t = time.time()
- nsteps = int(np.ceil((t_end + 1e-9) / dt))
- for (_, _, profile_data) in stepper.run( # noqa: B007
- ic, t_start, dt, t_end, return_profile_data=True):
- step += 1
- tn = time.time()
- logger.info("step %d/%d: %f", step, nsteps, tn - t)
- t = tn
-
- logger.info("fusion? %s", use_fusion)
- for key, value in profile_data.items():
- if isinstance(value, TimingFutureList):
- print(key, value.elapsed())
-
-# }}}
-
-
-# {{{ paper outputs
-
-def get_example_stepper(actx, dims=2, order=3, use_fusion=True,
- exec_mapper_factory=ExecutionMapper,
- return_ic=False):
- sym_operator, discr = get_wave_op_with_discr(
- actx, dims=dims, order=3)
-
- if not use_fusion:
- bound_op = bind(
- discr, sym_operator,
- exec_mapper_factory=exec_mapper_factory)
-
- stepper = RK4TimeStepper(
- discr, "w", bound_op, 1 + discr.dim,
- get_wave_component,
- exec_mapper_factory=exec_mapper_factory)
-
- else:
- stepper = FusedRK4TimeStepper(
- discr, "w", sym_operator, 1 + discr.dim,
- get_wave_component,
- exec_mapper_factory=exec_mapper_factory)
-
- if return_ic:
- ic = flat_obj_array(discr.zeros(actx),
- [discr.zeros(actx) for i in range(discr.dim)])
- return stepper, ic
-
- return stepper
-
-
-def latex_table(table_format, header, rows):
- result = []
- _ = result.append
- _(rf"\begin{{tabular}}{{{table_format}}}")
- _(r"\toprule")
- _(" & ".join(rf"\multicolumn{{1}}{{c}}{{{item}}}" for item in header) + r" \\")
- _(r"\midrule")
- for row in rows:
- _(" & ".join(row) + r" \\")
- _(r"\bottomrule")
- _(r"\end{tabular}")
- return "\n".join(result)
-
-
-def ascii_table(table_format, header, rows):
- from pytools import Table
- table = Table()
- table.add_row(header)
-
- for input_row in rows:
- row = []
- for item in input_row:
- if item.startswith(r"\num{"):
- # Strip \num{...} formatting
- row.append(item[5:-1])
- else:
- row.append(item)
- table.add_row(row)
-
- return str(table)
-
-
-if not PAPER_OUTPUT:
- table = ascii_table
-else:
- table = latex_table
-
-
-def problem_stats(order=3):
- cl_ctx = cl.create_some_context()
- queue = cl.CommandQueue(cl_ctx)
- actx = PyOpenCLArrayContext(
- queue,
- allocator=cl_tools.MemoryPool(cl_tools.ImmediateAllocator(queue))
- )
-
- with open_output_file("grudge-problem-stats.txt") as outf:
- _, dg_discr_2d = get_wave_op_with_discr(
- actx, dims=2, order=order)
- print("Number of 2D elements:", dg_discr_2d.mesh.nelements, file=outf)
- vol_discr_2d = dg_discr_2d.discr_from_dd("vol")
- dofs_2d = {group.nunit_dofs for group in vol_discr_2d.groups}
- from pytools import one
- print("Number of DOFs per 2D element:", one(dofs_2d), file=outf)
-
- _, dg_discr_3d = get_wave_op_with_discr(
- actx, dims=3, order=order)
- print("Number of 3D elements:", dg_discr_3d.mesh.nelements, file=outf)
- vol_discr_3d = dg_discr_3d.discr_from_dd("vol")
- dofs_3d = {group.nunit_dofs for group in vol_discr_3d.groups}
- from pytools import one
- print("Number of DOFs per 3D element:", one(dofs_3d), file=outf)
-
- logger.info("Wrote '%s'", outf.name)
-
-
-def statement_counts_table():
- cl_ctx = cl.create_some_context()
- queue = cl.CommandQueue(cl_ctx)
- actx = PyOpenCLArrayContext(
- queue,
- allocator=cl_tools.MemoryPool(cl_tools.ImmediateAllocator(queue))
- )
-
- fused_stepper = get_example_stepper(actx, use_fusion=True)
- stepper = get_example_stepper(actx, use_fusion=False)
-
- with open_output_file("statement-counts.tex") as outf:
- if not PAPER_OUTPUT:
- print("==== Statement Counts ====", file=outf)
-
- print(table(
- "lr",
- ("Operator", "Grudge Node Count"),
- (
- ("Time integration: baseline",
- r"\num{%d}"
- % len(stepper.bound_op.eval_code.instructions)),
- ("Right-hand side: baseline",
- r"\num{%d}"
- % len(stepper.grudge_bound_op.eval_code.instructions)),
- ("Inlined operator",
- r"\num{%d}"
- % len(fused_stepper.bound_op.eval_code.instructions))
- )),
- file=outf)
-
- logger.info("Wrote '%s'", outf.name)
-
-
-@memoize(key=lambda queue, dims: dims)
-def mem_ops_results(actx, dims):
- fused_stepper = get_example_stepper(
- actx,
- dims=dims,
- use_fusion=True,
- exec_mapper_factory=ExecutionMapperWithMemOpCounting)
-
- stepper, ic = get_example_stepper(
- actx,
- dims=dims,
- use_fusion=False,
- exec_mapper_factory=ExecutionMapperWithMemOpCounting,
- return_ic=True)
-
- t_start = 0
- dt = 0.02
- t_end = 0.02
-
- result = {}
-
- for (_, _, profile_data) in stepper.run( # noqa: B007
- ic, t_start, dt, t_end, return_profile_data=True):
- pass
-
- result["nonfused_bytes_read"] = profile_data["bytes_read"]
- result["nonfused_bytes_written"] = profile_data["bytes_written"]
- result["nonfused_bytes_total"] = \
- result["nonfused_bytes_read"] \
- + result["nonfused_bytes_written"]
-
- result["nonfused_bytes_read_by_scalar_assignments"] = \
- profile_data["bytes_read_by_scalar_assignments"]
- result["nonfused_bytes_written_by_scalar_assignments"] = \
- profile_data["bytes_written_by_scalar_assignments"]
- result["nonfused_bytes_total_by_scalar_assignments"] = \
- result["nonfused_bytes_read_by_scalar_assignments"] \
- + result["nonfused_bytes_written_by_scalar_assignments"]
-
- for (_, _, profile_data) in fused_stepper.run( # noqa: B007
- ic, t_start, dt, t_end, return_profile_data=True):
- pass
-
- result["fused_bytes_read"] = profile_data["bytes_read"]
- result["fused_bytes_written"] = profile_data["bytes_written"]
- result["fused_bytes_total"] = \
- result["fused_bytes_read"] \
- + result["fused_bytes_written"]
-
- result["fused_bytes_read_by_scalar_assignments"] = \
- profile_data["bytes_read_by_scalar_assignments"]
- result["fused_bytes_written_by_scalar_assignments"] = \
- profile_data["bytes_written_by_scalar_assignments"]
- result["fused_bytes_total_by_scalar_assignments"] = \
- result["fused_bytes_read_by_scalar_assignments"] \
- + result["fused_bytes_written_by_scalar_assignments"]
-
- return result
-
-
-def scalar_assignment_percent_of_total_mem_ops_table():
- cl_ctx = cl.create_some_context()
- queue = cl.CommandQueue(cl_ctx)
- actx = PyOpenCLArrayContext(
- queue,
- allocator=cl_tools.MemoryPool(cl_tools.ImmediateAllocator(queue))
- )
-
- result2d = mem_ops_results(actx, 2)
- result3d = mem_ops_results(actx, 3)
-
- with open_output_file("scalar-assignments-mem-op-percentage.tex") as outf:
- if not PAPER_OUTPUT:
- print("==== Scalar Assigment % of Total Mem Ops ====", file=outf)
-
- print(
- table(
- "lr",
- ("Operator",
- r"\parbox{1in}{\centering \% Memory Ops. "
- r"Due to Scalar Assignments}"),
- (
- ("2D: Baseline",
- "%.1f" % (
- 100 * result2d["nonfused_bytes_total_by_scalar_assignments"]
- / result2d["nonfused_bytes_total"])),
- ("2D: Inlined",
- "%.1f" % (
- 100 * result2d["fused_bytes_total_by_scalar_assignments"]
- / result2d["fused_bytes_total"])),
- ("3D: Baseline",
- "%.1f" % (
- 100 * result3d["nonfused_bytes_total_by_scalar_assignments"]
- / result3d["nonfused_bytes_total"])),
- ("3D: Inlined",
- "%.1f" % (
- 100 * result3d["fused_bytes_total_by_scalar_assignments"]
- / result3d["fused_bytes_total"])),
- )),
- file=outf)
-
- logger.info("Wrote '%s'", outf.name)
-
-
-def scalar_assignment_effect_of_fusion_mem_ops_table():
- cl_ctx = cl.create_some_context()
- queue = cl.CommandQueue(cl_ctx)
-
- result2d = mem_ops_results(queue, 2)
- result3d = mem_ops_results(queue, 3)
-
- with open_output_file("scalar-assignments-fusion-impact.tex") as outf:
- if not PAPER_OUTPUT:
- print("==== Scalar Assigment Inlining Impact ====", file=outf)
-
- print(
- table(
- "lrrrr",
- ("Operator",
- r"Bytes Read",
- r"Bytes Written",
- r"Total",
- r"\% of Baseline"),
- (
- ("2D: Baseline",
- r"\num{%d}" % (
- result2d["nonfused_bytes_read_by_scalar_assignments"]),
- r"\num{%d}" % (
- result2d["nonfused_bytes_written_by_scalar_assignments"]),
- r"\num{%d}" % (
- result2d["nonfused_bytes_total_by_scalar_assignments"]),
- "100"),
- ("2D: Inlined",
- r"\num{%d}" % (
- result2d["fused_bytes_read_by_scalar_assignments"]),
- r"\num{%d}" % (
- result2d["fused_bytes_written_by_scalar_assignments"]),
- r"\num{%d}" % (
- result2d["fused_bytes_total_by_scalar_assignments"]),
- r"%.1f" % (
- 100 * result2d["fused_bytes_total_by_scalar_assignments"]
- / result2d["nonfused_bytes_total_by_scalar_assignments"])),
- ("3D: Baseline",
- r"\num{%d}" % (
- result3d["nonfused_bytes_read_by_scalar_assignments"]),
- r"\num{%d}" % (
- result3d["nonfused_bytes_written_by_scalar_assignments"]),
- r"\num{%d}" % (
- result3d["nonfused_bytes_total_by_scalar_assignments"]),
- "100"),
- ("3D: Inlined",
- r"\num{%d}" % (
- result3d["fused_bytes_read_by_scalar_assignments"]),
- r"\num{%d}" % (
- result3d["fused_bytes_written_by_scalar_assignments"]),
- r"\num{%d}" % (
- result3d["fused_bytes_total_by_scalar_assignments"]),
- r"%.1f" % (
- 100 * result3d["fused_bytes_total_by_scalar_assignments"]
- / result3d["nonfused_bytes_total_by_scalar_assignments"])),
- )),
- file=outf)
- logger.info("Wrote '%s'", outf.name)
-
-# }}}
-
-
-def main():
- import sys
- if len(sys.argv) > 1:
- exec(sys.argv[1])
- else:
- if not SKIP_TESTS:
- # Run tests.
- result = pytest.main([__file__])
- assert result == 0
-
- # Run examples.
- problem_stats()
- statement_counts_table()
- scalar_assignment_percent_of_total_mem_ops_table()
- scalar_assignment_effect_of_fusion_mem_ops_table()
-
-
-if __name__ == "__main__":
- main()
-
-# vim: foldmethod=marker
diff --git a/examples/old_symbolics/symbolic_wave_op.py b/examples/old_symbolics/symbolic_wave_op.py
deleted file mode 100644
index 466c92a3..00000000
--- a/examples/old_symbolics/symbolic_wave_op.py
+++ /dev/null
@@ -1,166 +0,0 @@
-"""(Old world) symbolic wave equation operator."""
-
-__copyright__ = "Copyright (C) 2009 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
-from meshmode.mesh import BTAG_ALL, BTAG_NONE
-from grudge import sym
-from pytools.obj_array import flat_obj_array
-
-
-# {{{ Old-world symbolic wave operator
-
-class WeakWaveOperator:
- r"""This operator discretizes the wave equation
- :math:`\partial_t^2 u = c^2 \Delta u`.
-
- To be precise, we discretize the hyperbolic system
-
- .. math::
-
- \partial_t u - c \\nabla \\cdot v = 0
-
- \partial_t v - c \\nabla u = 0
-
- The sign of :math:`v` determines whether we discretize the forward or the
- backward wave equation.
-
- :math:`c` is assumed to be constant across all space.
- """
-
- def __init__(self, c, ambient_dim, source_f=0,
- flux_type="upwind",
- dirichlet_tag=BTAG_ALL,
- dirichlet_bc_f=0,
- neumann_tag=BTAG_NONE,
- radiation_tag=BTAG_NONE):
- assert isinstance(ambient_dim, int)
-
- self.c = c
- self.ambient_dim = ambient_dim
- self.source_f = source_f
-
- if self.c > 0:
- self.sign = 1
- else:
- self.sign = -1
-
- self.dirichlet_tag = dirichlet_tag
- self.neumann_tag = neumann_tag
- self.radiation_tag = radiation_tag
-
- self.dirichlet_bc_f = dirichlet_bc_f
-
- self.flux_type = flux_type
-
- def flux(self, w):
- u = w[0]
- v = w[1:]
- normal = sym.normal(w.dd, self.ambient_dim)
-
- central_flux_weak = -self.c*flat_obj_array(
- np.dot(v.avg, normal),
- u.avg * normal)
-
- if self.flux_type == "central":
- return central_flux_weak
- elif self.flux_type == "upwind":
- return central_flux_weak - self.c*self.sign*flat_obj_array(
- 0.5*(u.ext-u.int),
- 0.5*(normal * np.dot(normal, v.ext-v.int)))
- else:
- raise ValueError("invalid flux type '%s'" % self.flux_type)
-
- def sym_operator(self):
- d = self.ambient_dim
-
- w = sym.make_sym_array("w", d+1)
- u = w[0]
- v = w[1:]
-
- # boundary conditions -------------------------------------------------
-
- # dirichlet BCs -------------------------------------------------------
- dir_u = sym.cse(sym.project("vol", self.dirichlet_tag)(u))
- dir_v = sym.cse(sym.project("vol", self.dirichlet_tag)(v))
- if self.dirichlet_bc_f:
- # FIXME
- from warnings import warn
- warn("Inhomogeneous Dirichlet conditions on the wave equation "
- "are still having issues.")
-
- dir_g = sym.var("dir_bc_u")
- dir_bc = flat_obj_array(2*dir_g - dir_u, dir_v)
- else:
- dir_bc = flat_obj_array(-dir_u, dir_v)
-
- dir_bc = sym.cse(dir_bc, "dir_bc")
-
- # neumann BCs ---------------------------------------------------------
- neu_u = sym.cse(sym.project("vol", self.neumann_tag)(u))
- neu_v = sym.cse(sym.project("vol", self.neumann_tag)(v))
- neu_bc = sym.cse(flat_obj_array(neu_u, -neu_v), "neu_bc")
-
- # radiation BCs -------------------------------------------------------
- rad_normal = sym.normal(self.radiation_tag, d)
-
- rad_u = sym.cse(sym.project("vol", self.radiation_tag)(u))
- rad_v = sym.cse(sym.project("vol", self.radiation_tag)(v))
-
- rad_bc = sym.cse(flat_obj_array(
- 0.5*(rad_u - self.sign*np.dot(rad_normal, rad_v)),
- 0.5*rad_normal*(np.dot(rad_normal, rad_v) - self.sign*rad_u)
- ), "rad_bc")
-
- # entire operator -----------------------------------------------------
- def flux(pair):
- return sym.project(pair.dd, "all_faces")(self.flux(pair))
-
- result = sym.InverseMassOperator()(
- flat_obj_array(
- -self.c*np.dot(sym.stiffness_t(self.ambient_dim), v),
- -self.c*(sym.stiffness_t(self.ambient_dim)*u)
- )
-
- - sym.FaceMassOperator()(flux(sym.int_tpair(w))
- + flux(sym.bv_tpair(self.dirichlet_tag, w, dir_bc))
- + flux(sym.bv_tpair(self.neumann_tag, w, neu_bc))
- + flux(sym.bv_tpair(self.radiation_tag, w, rad_bc))
-
- ))
-
- result[0] += self.source_f
-
- return result
-
- 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])
-
-# }}}
-
-
-# vim: foldmethod=marker
diff --git a/grudge/__init__.py b/grudge/__init__.py
index 04cccfb1..aad8dbd1 100644
--- a/grudge/__init__.py
+++ b/grudge/__init__.py
@@ -20,11 +20,8 @@ OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN
THE SOFTWARE.
"""
-import grudge.symbolic as sym
-from grudge.execution import bind
-
from grudge.discretization import DiscretizationCollection
__all__ = [
- "sym", "bind", "DiscretizationCollection"
+ "DiscretizationCollection"
]
diff --git a/grudge/discretization.py b/grudge/discretization.py
index 017d91d7..3a86b706 100644
--- a/grudge/discretization.py
+++ b/grudge/discretization.py
@@ -158,17 +158,6 @@ class DiscretizationCollection:
self.group_factory_for_discretization_tag(DISCR_TAG_BASE)
)
- # NOTE: Can be removed when symbolics are completely removed
- # {{{ management of discretization-scoped common subexpressions
-
- from pytools import UniqueNameGenerator
- self._discr_scoped_name_gen = UniqueNameGenerator()
-
- self._discr_scoped_subexpr_to_name = {}
- self._discr_scoped_subexpr_name_to_value = {}
-
- # }}}
-
# {{{ process mpi_communicator argument
if mpi_communicator is not None:
diff --git a/grudge/execution.py b/grudge/execution.py
deleted file mode 100644
index 73b5e894..00000000
--- a/grudge/execution.py
+++ /dev/null
@@ -1,822 +0,0 @@
-__copyright__ = "Copyright (C) 2015-2017 Andreas Kloeckner, Bogdan Enache"
-
-__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 arraycontext import ArrayContext, make_loopy_program, thaw
-
-from typing import Optional, Union, Dict
-from numbers import Number
-import numpy as np
-
-from pytools import memoize_in
-
-import loopy as lp
-import pyopencl.array # noqa
-
-from meshmode.dof_array import DOFArray, flatten, unflatten
-
-import grudge.symbolic.mappers as mappers
-from grudge import sym
-from grudge.function_registry import base_function_registry
-
-import logging
-logger = logging.getLogger(__name__)
-
-from loopy.version import LOOPY_USE_LANGUAGE_VERSION_2018_2 # noqa: F401
-
-
-ResultType = Union[DOFArray, Number]
-
-
-# {{{ exec mapper
-
-class ExecutionMapper(mappers.Evaluator,
- mappers.BoundOpMapperMixin,
- mappers.LocalOpReducerMixin):
- def __init__(self, array_context, context, bound_op):
- super().__init__(context)
- self.dcoll = bound_op.dcoll
- self.bound_op = bound_op
- self.function_registry = bound_op.function_registry
- self.array_context = array_context
-
- # {{{ expression mappings
-
- def map_ones(self, expr):
- if expr.dd.is_scalar():
- return 1
-
- discr = self.dcoll.discr_from_dd(expr.dd)
-
- result = discr.empty(self.array_context)
- for grp_ary in result:
- grp_ary.fill(1.0)
- return result
-
- def map_node_coordinate_component(self, expr):
- discr = self.dcoll.discr_from_dd(expr.dd)
- return thaw(
- discr.nodes(
- # only save volume nodes or boundary nodes
- # (but not nodes for interior face discretizations, which
- # are likely only used once to compute the normals)
- cached=(
- discr.ambient_dim == discr.dim
- or expr.dd.is_boundary_or_partition_interface()
- )
- )[expr.axis],
- self.array_context
- )
-
- def map_grudge_variable(self, expr):
- from numbers import Number
-
- value = self.context[expr.name]
- if not expr.dd.is_scalar() and isinstance(value, Number):
- discr = self.dcoll.discr_from_dd(expr.dd)
- ary = discr.empty(self.array_context)
- for grp_ary in ary:
- grp_ary.fill(value)
- value = ary
-
- return value
-
- def map_subscript(self, expr):
- value = super().map_subscript(expr)
-
- if isinstance(expr.aggregate, sym.Variable):
- dd = expr.aggregate.dd
-
- from numbers import Number
- if not dd.is_scalar() and isinstance(value, Number):
- discr = self.dcoll.discr_from_dd(dd)
- ary = discr.empty(self.array_context)
- for grp_ary in ary:
- grp_ary.fill(value)
- value = ary
- return value
-
- def map_call(self, expr):
- args = [self.rec(p) for p in expr.parameters]
- return self.function_registry[expr.function.name](self.array_context, *args)
-
- # }}}
-
- # {{{ elementwise reductions
-
- def _map_elementwise_reduction(self, op_name, field_expr, dd):
- @memoize_in(self.array_context,
- (ExecutionMapper, "elementwise_%s_prg" % op_name))
- def prg():
- return make_loopy_program(
- "{[iel, idof, jdof]: 0<=iel<nelements and 0<=idof, jdof<ndofs}",
- """
- result[iel, idof] = %s(jdof, operand[iel, jdof])
- """ % op_name,
- name="grudge_elementwise_%s" % op_name)
-
- field = self.rec(field_expr)
- discr = self.dcoll.discr_from_dd(dd)
- assert field.shape == (len(discr.groups),)
-
- result = discr.empty(self.array_context, dtype=field.entry_dtype)
- for grp in discr.groups:
- assert field[grp.index].shape == (grp.nelements, grp.nunit_dofs)
- self.array_context.call_loopy(
- prg(),
- operand=field[grp.index],
- result=result[grp.index])
-
- return result
-
- def map_elementwise_sum(self, op, field_expr):
- return self._map_elementwise_reduction("sum", field_expr, op.dd_in)
-
- def map_elementwise_min(self, op, field_expr):
- return self._map_elementwise_reduction("min", field_expr, op.dd_in)
-
- def map_elementwise_max(self, op, field_expr):
- return self._map_elementwise_reduction("max", field_expr, op.dd_in)
-
- # }}}
-
- # {{{ nodal reductions
-
- def map_nodal_sum(self, op, field_expr):
- actx = self.array_context
- return sum([actx.np.sum(grp_ary)
- for grp_ary in self.rec(field_expr)])
-
- def map_nodal_max(self, op, field_expr):
- from functools import reduce
- actx = self.array_context
- return reduce(lambda acc, grp_ary: actx.np.maximum(acc,
- actx.np.max(grp_ary)),
- self.rec(field_expr), -np.inf)
-
- def map_nodal_min(self, op, field_expr):
- from functools import reduce
- actx = self.array_context
- return reduce(lambda acc, grp_ary: actx.np.minimum(acc,
- actx.np.min(grp_ary)),
- self.rec(field_expr), np.inf)
-
- # }}}
-
- def map_if(self, expr):
- bool_crit = self.rec(expr.condition)
-
- if isinstance(bool_crit, DOFArray):
- # continues below
- pass
- elif isinstance(bool_crit, (bool, np.bool_, np.number)):
- if bool_crit:
- return self.rec(expr.then)
- else:
- return self.rec(expr.else_)
- else:
- raise TypeError(
- "Expected criterion to be of type np.number or DOFArray")
-
- assert isinstance(bool_crit, DOFArray)
- ngroups = len(bool_crit)
-
- from pymbolic import var
- iel = var("iel")
- idof = var("idof")
- subscript = (iel, idof)
-
- then = self.rec(expr.then)
- else_ = self.rec(expr.else_)
-
- import pymbolic.primitives as p
- var = p.Variable
-
- if isinstance(then, DOFArray):
- sym_then = var("a")[subscript]
-
- def get_then(igrp):
- return then[igrp]
- elif isinstance(then, np.number):
- sym_then = var("a")
-
- def get_then(igrp):
- return then
- else:
- raise TypeError(
- "Expected 'then' to be of type np.number or DOFArray")
-
- if isinstance(else_, DOFArray):
- sym_else = var("b")[subscript]
-
- def get_else(igrp):
- return else_[igrp]
- elif isinstance(else_, np.number):
- sym_else = var("b")
-
- def get_else(igrp):
- return else_
- else:
- raise TypeError(
- "Expected 'else' to be of type np.number or DOFArray")
-
- @memoize_in(self.array_context, (ExecutionMapper, "map_if_knl"))
- def knl(sym_then, sym_else):
- return make_loopy_program(
- "{[iel, idof]: 0<=iel<nelements and 0<=idof<nunit_dofs}",
- [
- lp.Assignment(var("out")[iel, idof],
- p.If(var("crit")[iel, idof], sym_then, sym_else))
- ])
-
- return DOFArray(self.array_context, tuple(
- self.array_context.call_loopy(
- knl(sym_then, sym_else),
- crit=bool_crit[igrp],
- a=get_then(igrp),
- b=get_else(igrp))
- for igrp in range(ngroups)))
-
- # {{{ elementwise linear operators
-
- def map_ref_diff_base(self, op, field_expr):
- raise NotImplementedError(
- "differentiation should be happening in batched form")
-
- def _elwise_linear_loopy_prg(self):
- @memoize_in(self.array_context, (ExecutionMapper, "elwise_linear_knl"))
- def prg():
- result = make_loopy_program(
- """{[iel, idof, j]:
- 0<=iel<nelements and
- 0<=idof<ndiscr_nodes_out and
- 0<=j<ndiscr_nodes_in}""",
- "result[iel, idof] = sum(j, mat[idof, j] * vec[iel, j])",
- name="elwise_linear")
-
- result = lp.tag_array_axes(result, "mat", "stride:auto,stride:auto")
- return result
-
- return prg()
-
- 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
-
- in_discr = self.dcoll.discr_from_dd(op.dd_in)
- out_discr = self.dcoll.discr_from_dd(op.dd_out)
-
- result = out_discr.empty(self.array_context, dtype=field.entry_dtype)
-
- prg = self._elwise_linear_loopy_prg()
- for in_grp, out_grp in zip(in_discr.groups, out_discr.groups):
- cache_key = "elwise_linear", in_grp, out_grp, op, field.entry_dtype
- try:
- matrix = self.bound_op.operator_data_cache[cache_key]
- except KeyError:
- matrix = self.array_context.freeze(
- self.array_context.from_numpy(
- np.asarray(
- op.matrix(out_grp, in_grp),
- dtype=field.entry_dtype)))
-
- self.bound_op.operator_data_cache[cache_key] = matrix
-
- self.array_context.call_loopy(
- prg,
- mat=matrix,
- result=result[out_grp.index],
- vec=field[in_grp.index])
-
- return result
-
- def map_projection(self, op, field_expr):
- conn = self.dcoll.connection_from_dds(op.dd_in, op.dd_out)
- return conn(self.rec(field_expr))
-
- def map_opposite_partition_face_swap(self, op, field_expr):
- assert op.dd_in == op.dd_out
- bdry_conn = self.dcoll.distributed_boundary_swap_connection(op.dd_in)
- remote_bdry_vec = self.rec(field_expr) # swapped by RankDataSwapAssign
- return bdry_conn(remote_bdry_vec)
-
- def map_opposite_interior_face_swap(self, op, field_expr):
- return self.dcoll.opposite_face_connection()(self.rec(field_expr))
-
- # }}}
-
- # {{{ face mass operator
-
- def map_ref_face_mass_operator(self, op, field_expr):
- field = self.rec(field_expr)
-
- from grudge.tools import is_zero
- if is_zero(field):
- return 0
-
- @memoize_in(self.array_context, (ExecutionMapper, "face_mass_knl"))
- def prg():
- return make_loopy_program(
- """{[iel,idof,f,j]:
- 0<=iel<nelements and
- 0<=f<nfaces and
- 0<=idof<nvol_nodes and
- 0<=j<nface_nodes}""",
- """
- result[iel,idof] = sum(f, sum(j, mat[idof, f, j] * vec[f, iel, j]))
- """,
- name="face_mass")
-
- all_faces_conn = self.dcoll.connection_from_dds("vol", op.dd_in)
- all_faces_discr = all_faces_conn.to_discr
- vol_discr = all_faces_conn.from_discr
-
- result = vol_discr.empty(self.array_context, dtype=field.entry_dtype)
-
- assert len(all_faces_discr.groups) == len(vol_discr.groups)
-
- for afgrp, volgrp in zip(all_faces_discr.groups, vol_discr.groups):
- cache_key = "face_mass", afgrp, op, field.entry_dtype
-
- nfaces = volgrp.mesh_el_group.nfaces
-
- try:
- matrix = self.bound_op.operator_data_cache[cache_key]
- except KeyError:
- matrix = op.matrix(afgrp, volgrp, field.entry_dtype)
- matrix = self.array_context.freeze(
- self.array_context.from_numpy(matrix))
-
- self.bound_op.operator_data_cache[cache_key] = matrix
-
- input_view = field[afgrp.index].reshape(
- nfaces, volgrp.nelements, afgrp.nunit_dofs)
- self.array_context.call_loopy(
- prg(),
- mat=matrix,
- result=result[volgrp.index],
- vec=input_view)
-
- return result
-
- def map_signed_face_ones(self, expr):
- from grudge.dof_desc import DOFDesc, DD_VOLUME
-
- assert expr.dd.is_trace()
- face_discr = self.dcoll.discr_from_dd(expr.dd)
- assert face_discr.dim == 0
-
- # NOTE: ignore quadrature_tags on expr.dd, since we only care about
- # the face_id here
- all_faces_conn = self.dcoll.connection_from_dds(
- DD_VOLUME,
- DOFDesc(expr.dd.domain_tag))
-
- field = face_discr.empty(self.array_context, dtype=self.dcoll.real_dtype)
- for grp_ary in field:
- grp_ary.fill(1)
-
- for igrp, grp in enumerate(all_faces_conn.groups):
- for batch in grp.batches:
- i = self.array_context.thaw(batch.to_element_indices)
- grp_field = field[igrp].reshape(-1)
- grp_field[i] = \
- (2.0 * (batch.to_element_face % 2) - 1.0) * grp_field[i]
-
- return field
-
- # }}}
-
- # {{{ instruction execution functions
-
- def map_insn_rank_data_swap(self, insn, profile_data=None):
- local_data = self.array_context.to_numpy(flatten(self.rec(insn.field)))
- comm = self.dcoll.mpi_communicator
-
- # print("Sending data to rank %d with tag %d"
- # % (insn.i_remote_rank, insn.send_tag))
- send_req = comm.Isend(local_data, insn.i_remote_rank, tag=insn.send_tag)
-
- remote_data_host = np.empty_like(local_data)
- recv_req = comm.Irecv(remote_data_host, insn.i_remote_rank, insn.recv_tag)
-
- return [], [
- MPIRecvFuture(
- array_context=self.array_context,
- bdry_discr=self.dcoll.discr_from_dd(insn.dd_out),
- recv_req=recv_req,
- insn_name=insn.name,
- remote_data_host=remote_data_host),
- MPISendFuture(send_req)]
-
- def map_insn_loopy_kernel(self, insn, profile_data=None):
- kdescr = insn.kernel_descriptor
- discr = self.dcoll.discr_from_dd(kdescr.governing_dd)
-
- dof_array_kwargs = {}
- other_kwargs = {}
-
- for name, expr in kdescr.input_mappings.items():
- v = self.rec(expr)
- if isinstance(v, DOFArray):
- dof_array_kwargs[name] = v
- else:
- other_kwargs[name] = v
-
- for name in kdescr.scalar_args():
- v = other_kwargs[name]
- if isinstance(v, (int, float)):
- other_kwargs[name] = discr.real_dtype.type(v)
- elif isinstance(v, complex):
- other_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)))
-
- result = {}
- for grp in discr.groups:
- kwargs = other_kwargs.copy()
- kwargs["nelements"] = grp.nelements
- kwargs["nunit_dofs"] = grp.nunit_dofs
-
- for name, ary in dof_array_kwargs.items():
- kwargs[name] = ary[grp.index]
-
- knl_result = self.array_context.call_loopy(
- kdescr.loopy_kernel, **kwargs)
-
- for name, val in knl_result.items():
- result.setdefault(name, []).append(val)
-
- result = {
- name: DOFArray(self.array_context, tuple(val))
- for name, val in result.items()}
-
- return list(result.items()), []
-
- def map_insn_assign(self, insn, profile_data=None):
- return [(name, self.rec(expr))
- for name, expr in zip(insn.names, insn.exprs)], []
-
- def map_insn_assign_to_discr_scoped(self, insn, profile_data=None):
- assignments = []
- for name, expr in zip(insn.names, insn.exprs):
- value = self.rec(expr)
- self.dcoll._discr_scoped_subexpr_name_to_value[name] = value
- assignments.append((name, value))
-
- return assignments, []
-
- def map_insn_assign_from_discr_scoped(self, insn, profile_data=None):
- return [(insn.name,
- self.dcoll._discr_scoped_subexpr_name_to_value[insn.name])], []
-
- def map_insn_diff_batch_assign(self, insn, profile_data=None):
- field = self.rec(insn.field)
- repr_op = insn.operators[0]
-
- assert repr_op.dd_in.domain_tag == repr_op.dd_out.domain_tag
-
- in_discr = self.dcoll.discr_from_dd(repr_op.dd_in)
- out_discr = self.dcoll.discr_from_dd(repr_op.dd_out)
-
- prg = self._elwise_linear_loopy_prg()
-
- result = []
- for name, op in zip(insn.names, insn.operators):
- group_results = []
- for in_grp, out_grp in zip(in_discr.groups, out_discr.groups):
- assert in_grp.nelements == out_grp.nelements
-
- # Cache operator
- cache_key = "diff_batch", in_grp, out_grp, tuple(insn.operators),\
- field.entry_dtype
- try:
- matrices_dev = self.bound_op.operator_data_cache[cache_key]
- except KeyError:
- matrices_dev = [self.array_context.from_numpy(mat)
- for mat in repr_op.matrices(out_grp, in_grp)]
- self.bound_op.operator_data_cache[cache_key] = matrices_dev
-
- group_results.append(self.array_context.call_loopy(
- prg,
- mat=matrices_dev[op.rst_axis],
- vec=field[in_grp.index])["result"])
-
- result.append(
- (name, DOFArray(self.array_context, tuple(group_results))))
-
- return result, []
-
- # }}}
-
-# }}}
-
-
-# {{{ futures
-
-class MPIRecvFuture:
- def __init__(self, array_context, bdry_discr, recv_req, insn_name,
- remote_data_host):
- self.array_context = array_context
- self.bdry_discr = bdry_discr
- self.receive_request = recv_req
- self.insn_name = insn_name
- self.remote_data_host = remote_data_host
-
- def is_ready(self):
- return self.receive_request.Test()
-
- def __call__(self):
- self.receive_request.Wait()
- actx = self.array_context
- remote_data = unflatten(self.array_context, self.bdry_discr,
- actx.from_numpy(self.remote_data_host))
- return [(self.insn_name, remote_data)], []
-
-
-class MPISendFuture:
- def __init__(self, send_request):
- self.send_request = send_request
-
- def is_ready(self):
- return self.send_request.Test()
-
- def __call__(self):
- self.send_request.Wait()
- return [], []
-
-# }}}
-
-
-# {{{ bound operator
-
-class BoundOperator:
- def __init__(self, dcoll, discr_code, eval_code, debug_flags,
- function_registry, exec_mapper_factory):
- self.dcoll = dcoll
- self.discr_code = discr_code
- self.eval_code = eval_code
- self.operator_data_cache = {}
- self.debug_flags = debug_flags
- self.function_registry = function_registry
- self.exec_mapper_factory = exec_mapper_factory
-
- def __str__(self):
- sep = 75 * "=" + "\n"
- return (
- sep
- + "DISCRETIZATION-SCOPE CODE\n"
- + sep
- + str(self.discr_code) + "\n"
- + sep
- + "PER-EVALUATION CODE\n"
- + sep
- + str(self.eval_code))
-
- def __call__(self, array_context: Optional[ArrayContext] = None,
- *, profile_data=None, log_quantities=None, **context):
- """
- :arg array_context: only needs to be supplied if no instances of
- :class:`~meshmode.dof_array.DOFArray` with a
- :class:`~meshmode.array_context.ArrayContext`
- are supplied as part of *context*.
- """
-
- # {{{ figure array context
-
- array_contexts = []
- if array_context is not None:
- if not isinstance(array_context, ArrayContext):
- raise TypeError(
- "first positional argument (if supplied) must be "
- "an ArrayContext")
-
- array_contexts.append(array_context)
- del array_context
-
- def look_for_array_contexts(ary):
- if isinstance(ary, DOFArray):
- if ary.array_context is not None:
- array_contexts.append(ary.array_context)
- elif isinstance(ary, np.ndarray) and ary.dtype.char == "O":
- for idx in np.ndindex(ary.shape):
- look_for_array_contexts(ary[idx])
- else:
- pass
-
- for val in context.values():
- look_for_array_contexts(val)
-
- if array_contexts:
- from pytools import is_single_valued
- if not is_single_valued(array_contexts):
- raise ValueError("arguments do not agree on an array context")
-
- array_context = array_contexts[0]
- else:
- raise ValueError("no array context given or available from arguments")
-
- # }}}
-
- # {{{ dcoll-scope evaluation
-
- if any(
- (result_var.name not in
- self.dcoll._discr_scoped_subexpr_name_to_value)
- for result_var in self.discr_code.result):
- # need to do dcoll-scope evaluation
- dcoll_eval_context: Dict[str, ResultType] = {}
- self.discr_code.execute(
- self.exec_mapper_factory(
- array_context, dcoll_eval_context, self))
-
- # }}}
-
- return self.eval_code.execute(
- self.exec_mapper_factory(array_context, context, self),
- profile_data=profile_data,
- log_quantities=log_quantities)
-
-# }}}
-
-
-# {{{ process_sym_operator function
-
-def process_sym_operator(dcoll, sym_operator, post_bind_mapper=None, dumper=None,
- local_only=None):
- if local_only is None:
- local_only = False
-
- if dumper is None:
- def dumper(name, sym_operator):
- return
-
- orig_sym_operator = sym_operator
- import grudge.symbolic.mappers as mappers
-
- dumper("before-bind", sym_operator)
- sym_operator = mappers.OperatorBinder()(sym_operator)
-
- mappers.ErrorChecker(dcoll.mesh)(sym_operator)
-
- sym_operator = \
- mappers.OppositeInteriorFaceSwapUniqueIDAssigner()(sym_operator)
-
- if not local_only:
- # {{{ broadcast root rank's symn_operator
-
- # also make sure all ranks had same orig_sym_operator
-
- if dcoll.mpi_communicator is not None:
- (mgmt_rank_orig_sym_operator, mgmt_rank_sym_operator) = \
- dcoll.mpi_communicator.bcast(
- (orig_sym_operator, sym_operator),
- dcoll.get_management_rank_index())
-
- if not np.array_equal(mgmt_rank_orig_sym_operator, orig_sym_operator):
- raise ValueError("rank %d received a different symbolic "
- "operator to bind from rank %d"
- % (dcoll.mpi_communicator.Get_rank(),
- dcoll.get_management_rank_index()))
-
- sym_operator = mgmt_rank_sym_operator
-
- # }}}
-
- if post_bind_mapper is not None:
- dumper("before-postbind", sym_operator)
- sym_operator = post_bind_mapper(sym_operator)
-
- dumper("before-empty-flux-killer", sym_operator)
- sym_operator = mappers.EmptyFluxKiller(dcoll.mesh)(sym_operator)
-
- dumper("before-cfold", sym_operator)
- sym_operator = mappers.CommutativeConstantFoldingMapper()(sym_operator)
-
- dumper("before-qcheck", sym_operator)
- sym_operator = mappers.QuadratureCheckerAndRemover(
- dcoll.discr_tag_to_group_factory)(sym_operator)
-
- # Work around https://github.com/numpy/numpy/issues/9438
- #
- # The idea is that we need 1j as an expression to survive
- # until code generation time. If it is evaluated and combined
- # with other constants, we will need to determine its size
- # (as np.complex64/128) within the expression. But because
- # of the above numpy bug, sized numbers are not likely to survive
- # expression building--so that's why we step in here to fix that.
-
- dumper("before-csize", sym_operator)
- sym_operator = mappers.ConstantToNumpyConversionMapper(
- real_type=dcoll.real_dtype.type,
- complex_type=dcoll.complex_dtype.type,
- )(sym_operator)
-
- dumper("before-global-to-reference", sym_operator)
- sym_operator = mappers.GlobalToReferenceMapper(dcoll)(sym_operator)
-
- dumper("before-distributed", sym_operator)
-
- if not local_only:
- volume_mesh = dcoll.discr_from_dd("vol").mesh
- from meshmode.distributed import get_connected_partitions
- connected_parts = get_connected_partitions(volume_mesh)
-
- if connected_parts:
- sym_operator = mappers.DistributedMapper(connected_parts)(sym_operator)
-
- dumper("before-imass", sym_operator)
- sym_operator = mappers.InverseMassContractor()(sym_operator)
-
- dumper("before-cfold-2", sym_operator)
- sym_operator = mappers.CommutativeConstantFoldingMapper()(sym_operator)
-
- # FIXME: Reenable derivative joiner
- # dumper("before-derivative-join", sym_operator)
- # sym_operator = mappers.DerivativeJoiner()(sym_operator)
-
- dumper("process-finished", sym_operator)
-
- return sym_operator
-
-# }}}
-
-
-def bind(discr, sym_operator, *, post_bind_mapper=lambda x: x,
- function_registry=base_function_registry,
- exec_mapper_factory=ExecutionMapper,
- debug_flags=frozenset(), local_only=None):
- """
- :param local_only: If *True*, *sym_operator* should oly be evaluated on the
- local part of the mesh. No inter-rank communication will take place.
- (However rank boundaries, tagged :class:`~meshmode.mesh.BTAG_PARTITION`,
- will not automatically be considered part of the domain boundary.)
- """
- # from grudge.symbolic.mappers import QuadratureUpsamplerRemover
- # sym_operator = QuadratureUpsamplerRemover(self.quad_min_degrees)(
- # sym_operator)
-
- stage = [0]
-
- def dump_sym_operator(name, sym_operator):
- if "dump_sym_operator_stages" in debug_flags:
- from pytools.debug import open_unique_debug_file
- outf, name = open_unique_debug_file("%02d-%s" % (stage[0], name), ".txt")
- with outf:
- outf.write(sym.pretty(sym_operator))
-
- stage[0] += 1
-
- sym_operator = process_sym_operator(
- discr,
- sym_operator,
- post_bind_mapper=post_bind_mapper,
- dumper=dump_sym_operator,
- local_only=local_only)
-
- from grudge.symbolic.compiler import OperatorCompiler
- discr_code, eval_code = OperatorCompiler(discr, function_registry)(sym_operator)
-
- bound_op = BoundOperator(discr, discr_code, eval_code,
- function_registry=function_registry,
- exec_mapper_factory=exec_mapper_factory,
- debug_flags=debug_flags)
-
- if "dump_op_code" in debug_flags:
- from pytools.debug import open_unique_debug_file
- outf, _ = open_unique_debug_file("op-code", ".txt")
- with outf:
- outf.write(str(bound_op))
-
- if "dump_dataflow_graph" in debug_flags:
- discr_code.dump_dataflow_graph("discr")
- eval_code.dump_dataflow_graph("eval")
-
- return bound_op
-
-# vim: foldmethod=marker
diff --git a/grudge/function_registry.py b/grudge/function_registry.py
deleted file mode 100644
index 89629def..00000000
--- a/grudge/function_registry.py
+++ /dev/null
@@ -1,193 +0,0 @@
-__copyright__ = """
-Copyright (C) 2013 Andreas Kloeckner
-Copyright (C) 2019 Matt Wala
-"""
-
-__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
-
-from loopy.version import LOOPY_USE_LANGUAGE_VERSION_2018_2 # noqa
-from pytools import RecordWithoutPickling
-
-
-# {{{ helpers
-
-def should_use_numpy(arg):
- from numbers import Number
- if isinstance(arg, Number) or \
- isinstance(arg, np.ndarray) and arg.shape == ():
- return True
- return False
-
-
-def cl_to_numpy_function_name(name):
- return {
- "atan2": "arctan2",
- }.get(name, name)
-
-# }}}
-
-
-# {{{ function
-
-class FunctionNotFound(KeyError):
- pass
-
-
-class Function(RecordWithoutPickling):
- """
- .. attribute:: identifier
- .. attribute:: supports_codegen
- .. automethod:: __call__
- .. automethod:: get_result_dofdesc
- """
-
- def __init__(self, identifier, **kwargs):
- super().__init__(identifier=identifier, **kwargs)
-
- def __call__(self, queue, *args, **kwargs):
- """Call the function implementation, if available."""
- raise TypeError("function '%s' is not callable" % self.identifier)
-
- def get_result_dofdesc(self, arg_dds):
- """Return the :class:`grudge.symbolic.primitives.DOFDesc` for the return value
- of the function.
-
- :arg arg_dds: A list of :class:`grudge.symbolic.primitives.DOFDesc` instances
- for each argument
- """
- raise NotImplementedError
-
-
-class CElementwiseFunction(Function):
- supports_codegen = True
-
- def __init__(self, identifier, nargs):
- super().__init__(identifier=identifier, nargs=nargs)
-
- def get_result_dofdesc(self, arg_dds):
- assert len(arg_dds) == self.nargs
- return arg_dds[0]
-
- def __call__(self, array_context, *args):
- func_name = self.identifier
- from pytools import single_valued
- if single_valued(should_use_numpy(arg) for arg in args):
- func = getattr(np, func_name)
- return func(*args)
-
- if func_name == "fabs": # FIXME
- # Loopy has a type-adaptive "abs", but no "fabs".
- func_name = "abs"
-
- sfunc = getattr(array_context.np, func_name)
- return sfunc(*args)
-
-
-class CBesselFunction(Function):
-
- supports_codegen = True
-
- def get_result_dofdesc(self, arg_dds):
- assert len(arg_dds) == 2
- return arg_dds[1]
-
-
-class FixedDOFDescExternalFunction(Function):
-
- supports_codegen = False
-
- def __init__(self, identifier, implementation, dd):
- super().__init__(
- identifier,
- implementation=implementation,
- dd=dd)
-
- def __call__(self, array_context, *args, **kwargs):
- return self.implementation(array_context, *args, **kwargs)
-
- def get_result_dofdesc(self, arg_dds):
- return self.dd
-
-# }}}
-
-
-# {{{ function registry
-
-class FunctionRegistry(RecordWithoutPickling):
- def __init__(self, id_to_function=None):
- if id_to_function is None:
- id_to_function = {}
-
- super().__init__(id_to_function=id_to_function)
-
- def register(self, function):
- """Return a copy of *self* with *function* registered."""
-
- if function.identifier in self.id_to_function:
- raise ValueError("function '%s' is already registered"
- % function.identifier)
-
- new_id_to_function = self.id_to_function.copy()
- new_id_to_function[function.identifier] = function
- return self.copy(id_to_function=new_id_to_function)
-
- def __getitem__(self, function_id):
- try:
- return self.id_to_function[function_id]
- except KeyError:
- raise FunctionNotFound(
- "unknown function: '%s'"
- % function_id)
-
- def __contains__(self, function_id):
- return function_id in self.id_to_function
-
-# }}}
-
-
-def _make_bfr():
- bfr = FunctionRegistry()
-
- bfr = bfr.register(CElementwiseFunction("sqrt", 1))
- bfr = bfr.register(CElementwiseFunction("exp", 1))
- bfr = bfr.register(CElementwiseFunction("fabs", 1))
- bfr = bfr.register(CElementwiseFunction("sin", 1))
- bfr = bfr.register(CElementwiseFunction("cos", 1))
- bfr = bfr.register(CElementwiseFunction("atan2", 2))
- bfr = bfr.register(CBesselFunction("bessel_j"))
- bfr = bfr.register(CBesselFunction("bessel_y"))
-
- return bfr
-
-
-base_function_registry = _make_bfr()
-
-
-def register_external_function(
- function_registry, identifier, implementation, dd):
- return function_registry.register(
- FixedDOFDescExternalFunction(
- identifier, implementation, dd))
-
-# vim: foldmethod=marker
diff --git a/grudge/symbolic/__init__.py b/grudge/symbolic/__init__.py
deleted file mode 100644
index c0756740..00000000
--- a/grudge/symbolic/__init__.py
+++ /dev/null
@@ -1,29 +0,0 @@
-"""Building blocks and mappers for operator expression trees."""
-
-__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 grudge.symbolic.primitives import * # noqa: F401,F403
-from grudge.symbolic.operators import * # noqa: F401,F403
-from grudge.symbolic.tools import pretty # noqa: F401
-
-from pymbolic.primitives import If, Comparison # noqa: F401
diff --git a/grudge/symbolic/compiler.py b/grudge/symbolic/compiler.py
deleted file mode 100644
index 207ae189..00000000
--- a/grudge/symbolic/compiler.py
+++ /dev/null
@@ -1,1369 +0,0 @@
-"""Compiler to turn operator expression tree into (imperative) bytecode."""
-
-__copyright__ = "Copyright (C) 2008-15 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
-
-from pytools import Record, memoize_method, memoize
-from pytools.obj_array import obj_array_vectorize
-
-from grudge import sym
-import grudge.symbolic.mappers as mappers
-from pymbolic.primitives import Variable, Subscript
-from sys import intern
-from functools import reduce
-
-from loopy import ScalarCallable
-from loopy.version import LOOPY_USE_LANGUAGE_VERSION_2018_2 # noqa: F401
-
-
-# {{{ instructions
-
-class Instruction(Record):
- priority = 0
- neglect_for_dofdesc_inference = False
-
- def get_assignees(self):
- raise NotImplementedError("no get_assignees in %s" % self.__class__)
-
- def get_dependencies(self):
- raise NotImplementedError("no get_dependencies in %s" % self.__class__)
-
- def __str__(self):
- raise NotImplementedError
-
- def __hash__(self):
- return id(self)
-
- def __eq__(self, other):
- return self is other
-
- def __ne__(self, other):
- return not self.__eq__(other)
-
-
-@memoize(use_kwargs=True)
-def _make_dep_mapper(include_subscripts):
- return mappers.DependencyMapper(
- include_operator_bindings=False,
- include_subscripts=include_subscripts,
- include_calls="descend_args")
-
-
-# {{{ loopy kernel instruction
-
-class LoopyKernelDescriptor:
- 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.default_entrypoint.args
- if isinstance(arg, lp.ValueArg)
- and arg.name not in ["nelements", "nunit_dofs"]]
-
-
-class LoopyKernelInstruction(Instruction):
- comment = ""
- scope_indicator = ""
-
- def __init__(self, kernel_descriptor):
- super().__init__()
- self.kernel_descriptor = kernel_descriptor
-
- @memoize_method
- def get_assignees(self):
- return {k for k in self.kernel_descriptor.output_mappings.keys()}
-
- @memoize_method
- def get_dependencies(self):
- from pymbolic.primitives import Variable, Subscript
- result = set()
- for v in self.kernel_descriptor.input_mappings.values():
- if isinstance(v, Variable):
- result.add(v)
- elif isinstance(v, Subscript):
- result.add(v.aggregate)
- return result
-
- def __str__(self):
- knl_str = "\n".join(
- f"{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 = ""
-
- def __str__(self):
- comment = self.comment
- if len(self.names) == 1:
- if comment:
- comment = "/* %s */ " % comment
-
- return "{} <-{} {}{}".format(
- self.names[0], self.scope_indicator, comment,
- self.exprs[0])
- else:
- if comment:
- comment = " /* %s */" % comment
-
- lines = []
- lines.append("{" + comment)
- for n, e, dnr in zip(self.names, self.exprs, self.do_not_return):
- if dnr:
- dnr_indicator = "-#"
- else:
- dnr_indicator = ""
-
- lines.append(" {} <{}-{} {}".format(
- n, dnr_indicator, self.scope_indicator, e))
- lines.append("}")
- return "\n".join(lines)
-
-
-class Assign(AssignBase):
- """
- .. attribute:: names
- .. attribute:: exprs
- .. attribute:: do_not_return
-
- a list of bools indicating whether the corresponding entry in names and
- exprs describes an expression that is not needed beyond this assignment
-
- .. attribute:: priority
- """
-
- def __init__(self, names, exprs, **kwargs):
- Instruction.__init__(self, names=names, exprs=exprs, **kwargs)
-
- if not hasattr(self, "do_not_return"):
- self.do_not_return = [False] * len(names)
-
- @memoize_method
- def flop_count(self):
- return sum(mappers.FlopCounter()(expr) for expr in self.exprs)
-
- def get_assignees(self):
- return set(self.names)
-
- @memoize_method
- def get_dependencies(self, each_vector=False):
- dep_mapper = _make_dep_mapper(include_subscripts=False)
-
- from operator import or_
- deps = reduce(
- or_, (dep_mapper(expr)
- for expr in self.exprs))
-
- from pymbolic.primitives import Variable
- deps -= {Variable(name) for name in self.names}
-
- if not each_vector:
- self._dependencies = deps
-
- return deps
-
- mapper_method = intern("map_insn_assign")
-
-
-class RankDataSwapAssign(Instruction):
- """
- .. attribute:: name
- .. attribute:: field
- .. attribute:: i_remote_rank
-
- The number of the remote rank that this instruction swaps data with.
-
- .. attribute:: dd_out
- .. attribute:: comment
- """
- MPI_TAG_GRUDGE_DATA_BASE = 15165
-
- def __init__(self, name, field, op):
- self.name = name
- self.field = field
- self.i_remote_rank = op.i_remote_part
- self.dd_out = op.dd_out
- self.send_tag = self.MPI_TAG_GRUDGE_DATA_BASE + op.unique_id
- self.recv_tag = self.MPI_TAG_GRUDGE_DATA_BASE + op.unique_id
- self.comment = f"Swap data with rank {self.i_remote_rank:02d}"
-
- @memoize_method
- def get_assignees(self):
- return {self.name}
-
- @memoize_method
- def get_dependencies(self):
- return _make_dep_mapper(include_subscripts=False)(self.field)
-
- def __str__(self):
- return ("{\n"
- + f" /* {self.comment} */\n"
- + f" send_tag = {self.send_tag}\n"
- + f" recv_tag = {self.recv_tag}\n"
- + f" {self.name} <- {self.field}\n"
- + "}")
-
- mapper_method = intern("map_insn_rank_data_swap")
-
-
-class ToDiscretizationScopedAssign(Assign):
- scope_indicator = "(to discr)-"
-
- 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().__init__(name=name, **kwargs)
-
- @memoize_method
- def flop_count(self):
- return 0
-
- def get_assignees(self):
- return frozenset([self.name])
-
- def get_dependencies(self):
- return frozenset()
-
- def __str__(self):
- return "%s <-(from discr)" % self.name
-
- mapper_method = intern("map_insn_assign_from_discr_scoped")
-
-
-class DiffBatchAssign(Instruction):
- """
- :ivar names:
- :ivar operators:
-
- .. note ::
-
- All operators here are guaranteed to satisfy
- :meth:`grudge.symbolic.operators.DiffOperatorBase.
- equal_except_for_axis`.
-
- :ivar field:
- """
-
- def get_assignees(self):
- return frozenset(self.names)
-
- @memoize_method
- def get_dependencies(self):
- return _make_dep_mapper(include_subscripts=False)(self.field)
-
- def __str__(self):
- lines = []
-
- if len(self.names) > 1:
- lines.append("{")
- for n, d in zip(self.names, self.operators):
- lines.append(f" {n} <- {d}({self.field})")
- lines.append("}")
- else:
- for n, d in zip(self.names, self.operators):
- lines.append(f"{n} <- {d}({self.field})")
-
- return "\n".join(lines)
-
- mapper_method = intern("map_insn_diff_batch_assign")
-
-# }}}
-
-
-# {{{ graphviz/dot dataflow graph drawing
-
-def dot_dataflow_graph(code, max_node_label_length=30,
- label_wrap_width=50):
- origins = {}
- node_names = {}
-
- result = [
- 'initial [label="initial"]'
- 'result [label="result"]']
-
- for num, insn in enumerate(code.instructions):
- node_name = "node%d" % num
- node_names[insn] = node_name
- node_label = str(insn)
-
- 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:
- from pytools import word_wrap
- node_label = word_wrap(node_label, label_wrap_width,
- wrap_using="\n ")
-
- node_label = node_label.replace("\n", "\\l") + "\\l"
-
- result.append(f"{node_name} [ "
- f'label="p{insn.priority}: {node_label}" shape=box ];')
-
- for assignee in insn.get_assignees():
- origins[assignee] = node_name
-
- def get_orig_node(expr):
- from pymbolic.primitives import Variable
- if isinstance(expr, Variable):
- return origins.get(expr.name, "initial")
- else:
- return "initial"
-
- def gen_expr_arrow(expr, target_node):
- orig_node = get_orig_node(expr)
- result.append(f'{orig_node} -> {target_node} [label="{expr}"];')
-
- for insn in code.instructions:
- for dep in insn.get_dependencies():
- gen_expr_arrow(dep, node_names[insn])
-
- if isinstance(code.result, np.ndarray) and code.result.dtype.char == "O":
- for subexp in code.result:
- gen_expr_arrow(subexp, "result")
- else:
- gen_expr_arrow(code.result, "result")
-
- return "digraph dataflow {\n%s\n}\n" % "\n".join(result)
-
-# }}}
-
-
-# {{{ code representation
-
-class Code:
- def __init__(self, instructions, result):
- self.instructions = instructions
- self.result = result
- # self.last_schedule = None
- self.static_schedule_attempts = 5
-
- def dump_dataflow_graph(self, name=None):
- from pytools.debug import open_unique_debug_file
-
- if name is None:
- stem = "dataflow"
- else:
- stem = "dataflow-%s" % name
-
- outf, _ = open_unique_debug_file(stem, ".dot")
- with outf:
- outf.write(dot_dataflow_graph(self, max_node_label_length=None))
-
- def __str__(self):
- var_to_writer = {
- var_name: insn
- for insn in self.instructions
- for var_name in insn.get_assignees()}
-
- # {{{ topological sort
-
- added_insns = set()
- ordered_insns = []
-
- def insert_insn(insn):
- if insn in added_insns:
- return
-
- for dep in insn.get_dependencies():
- try:
- if isinstance(dep, Subscript):
- dep_name = dep.aggregate.name
- else:
- dep_name = dep.name
-
- writer = var_to_writer[dep_name]
- except KeyError:
- # input variables won't be found
- pass
- else:
- insert_insn(writer)
-
- ordered_insns.append(insn)
- added_insns.add(insn)
-
- for insn in self.instructions:
- insert_insn(insn)
-
- assert len(ordered_insns) == len(self.instructions)
- assert len(added_insns) == len(self.instructions)
-
- # }}}
-
- lines = []
- for insn in ordered_insns:
- lines.extend(str(insn).split("\n"))
- lines.append("RESULT: " + str(self.result))
-
- return "\n".join(lines)
-
- # {{{ dynamic scheduler (generates static schedules by self-observation)
-
- class NoInstructionAvailable(Exception):
- pass
-
- @memoize_method
- def get_next_step(self, available_names, done_insns):
- from pytools import argmax2
- available_insns = [
- (insn, insn.priority) for insn in self.instructions
- if insn not in done_insns
- and all(dep.name in available_names
- for dep in insn.get_dependencies())]
-
- if not available_insns:
- raise self.NoInstructionAvailable
-
- from pytools import flatten
- discardable_vars = set(available_names) - set(flatten(
- [dep.name for dep in insn.get_dependencies()]
- for insn in self.instructions
- if insn not in done_insns))
-
- # {{{ make sure results do not get discarded
-
- dm = mappers.DependencyMapper(composite_leaves=False)
-
- def remove_result_variable(result_expr):
- # The extra dependency mapper run is necessary
- # because, for instance, subscripts can make it
- # into the result expression, which then does
- # not consist of just variables.
-
- for var in dm(result_expr):
- assert isinstance(var, Variable)
- discardable_vars.discard(var.name)
-
- obj_array_vectorize(remove_result_variable, self.result)
-
- # }}}
-
- return argmax2(available_insns), discardable_vars
-
- def execute(self, exec_mapper, pre_assign_check=None, profile_data=None,
- log_quantities=None):
- if profile_data is not None:
- from time import time
- start_time = time()
- if profile_data == {}:
- profile_data["insn_eval_time"] = 0
- profile_data["future_eval_time"] = 0
- profile_data["busy_wait_time"] = 0
- profile_data["total_time"] = 0
- if log_quantities is not None:
- exec_sub_timer = log_quantities["exec_timer"].start_sub_timer()
- context = exec_mapper.context
-
- futures = []
- done_insns = set()
-
- while True:
- try:
- if profile_data is not None:
- insn_start_time = time()
- if log_quantities is not None:
- insn_sub_timer = \
- log_quantities["insn_eval_timer"].start_sub_timer()
-
- insn, discardable_vars = self.get_next_step(
- frozenset(list(context.keys())),
- frozenset(done_insns))
-
- done_insns.add(insn)
- for name in discardable_vars:
- del context[name]
-
- mapper_method = getattr(exec_mapper, insn.mapper_method)
- if log_quantities is not None:
- if isinstance(insn, RankDataSwapAssign):
- from logpyle import time_and_count_function
- mapper_method = time_and_count_function(
- mapper_method,
- log_quantities["rank_data_swap_timer"],
- log_quantities["rank_data_swap_counter"])
-
- assignments, new_futures = mapper_method(insn, profile_data)
-
- for target, value in assignments:
- if pre_assign_check is not None:
- pre_assign_check(target, value)
- context[target] = value
-
- futures.extend(new_futures)
- if profile_data is not None:
- profile_data["insn_eval_time"] += time() - insn_start_time
- if log_quantities is not None:
- insn_sub_timer.stop().submit()
- except self.NoInstructionAvailable:
- if not futures:
- # No more instructions or futures. We are done.
- break
-
- # Busy wait for a new future
- if profile_data is not None:
- busy_wait_start_time = time()
- if log_quantities is not None:
- busy_sub_timer =\
- log_quantities["busy_wait_timer"].start_sub_timer()
-
- did_eval_future = False
- while not did_eval_future:
- for i in range(len(futures)):
- if futures[i].is_ready():
- if profile_data is not None:
- profile_data["busy_wait_time"] +=\
- time() - busy_wait_start_time
- future_start_time = time()
- if log_quantities is not None:
- busy_sub_timer.stop().submit()
- future_sub_timer =\
- log_quantities["future_eval_timer"]\
- .start_sub_timer()
-
- future = futures.pop(i)
- assignments, new_futures = future()
-
- for target, value in assignments:
- if pre_assign_check is not None:
- pre_assign_check(target, value)
- context[target] = value
-
- futures.extend(new_futures)
- did_eval_future = True
-
- if profile_data is not None:
- profile_data["future_eval_time"] +=\
- time() - future_start_time
- if log_quantities is not None:
- future_sub_timer.stop().submit()
- break
-
- if len(done_insns) < len(self.instructions):
- raise RuntimeError("not all instructions are reachable"
- "--did you forget to pass a value for a placeholder?")
-
- if log_quantities is not None:
- exec_sub_timer.stop().submit()
- if profile_data is not None:
- profile_data["total_time"] = time() - start_time
- return (obj_array_vectorize(exec_mapper, self.result),
- profile_data)
- return obj_array_vectorize(exec_mapper, self.result)
-
-# }}}
-
-# }}}
-
-
-# {{{ assignment aggregration pass
-
-def aggregate_assignments(inf_mapper, instructions, result,
- max_vectors_in_batch_expr):
- from pymbolic.primitives import Variable
-
- function_registry = inf_mapper.function_registry
-
- # {{{ aggregation helpers
-
- def get_complete_origins_set(insn, skip_levels=0):
- try:
- return insn_to_origins_cache[insn]
- except KeyError:
- pass
-
- 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)
-
- insn_to_origins_cache[insn] = result
-
- return result
-
- var_assignees_cache = {}
-
- def get_var_assignees(insn):
- try:
- return var_assignees_cache[insn]
- except KeyError:
- result = {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()) \
- - {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
-
- insn_to_origins_cache = {}
-
- origins_map = {
- assignee: insn
- for insn in instructions
- for assignee in insn.get_assignees()}
-
- from pytools import partition
- from grudge.dof_desc 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 is_external_call(insn.exprs[0], function_registry)
- 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 = {
- expr
- for insn in processed_assigns + other_insns
- for expr in insn.get_dependencies()}
-
- if isinstance(result, np.ndarray) and result.dtype.char == "O":
- externally_used_names |= {expr for expr in result}
- else:
- externally_used_names |= {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 = {name for name, expr in names_exprs}
- names_exprs_deps = [
- (name, expr,
- {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
-
- # }}}
-
-# }}}
-
-
-# {{{ to-loopy mapper
-
-def is_external_call(expr, function_registry):
- from pymbolic.primitives import Call
- if not isinstance(expr, Call):
- return False
- return not is_function_loopyable(expr.function, function_registry)
-
-
-def is_function_loopyable(function, function_registry):
- from grudge.symbolic.primitives import FunctionSymbol
- assert isinstance(function, FunctionSymbol)
- return function_registry[function.name].supports_codegen
-
-
-class ToLoopyExpressionMapper(mappers.IdentityMapper):
- def __init__(self, dd_inference_mapper, temp_names, subscript):
- self.dd_inference_mapper = dd_inference_mapper
- self.function_registry = dd_inference_mapper.function_registry
- self.temp_names = temp_names
- self.subscript = subscript
-
- self.expr_to_name = {}
- self.used_names = set()
- self.non_scalar_vars = []
-
- def map_name(self, name):
- dot_idx = name.find(".")
- if dot_idx != -1:
- return "grdg_sub_{}_{}".format(name[:dot_idx], name[dot_idx+1:])
- else:
- return name
-
- def map_variable_ref_expr(self, expr, name_prefix):
- from pymbolic import var
- dd = self.dd_inference_mapper(expr)
-
- try:
- name = self.expr_to_name[expr]
- except KeyError:
- name_prefix = self.map_name(name_prefix)
- name = name_prefix
-
- suffix_nr = 0
- while name in self.used_names:
- name = f"{name_prefix}_{suffix_nr}"
- suffix_nr += 1
- self.used_names.add(name)
-
- self.expr_to_name[expr] = name
-
- from grudge.dof_desc import DTAG_SCALAR
- if dd.domain_tag == DTAG_SCALAR or name in self.temp_names:
- return var(name)
- else:
- self.non_scalar_vars.append(name)
- return var(name)[self.subscript]
-
- def map_variable(self, expr):
- return self.map_variable_ref_expr(expr, expr.name)
-
- def map_grudge_variable(self, expr):
- return self.map_variable_ref_expr(expr, expr.name)
-
- def map_subscript(self, expr):
- subscript = expr.index
- if isinstance(subscript, tuple):
- assert len(subscript) == 1
- subscript, = subscript
-
- assert isinstance(subscript, int)
-
- return self.map_variable_ref_expr(
- expr,
- "%s_%d" % (expr.aggregate.name, subscript))
-
- def map_call(self, expr):
- if is_function_loopyable(expr.function, self.function_registry):
- from pymbolic import var
-
- func_name = expr.function.name
- if func_name == "fabs":
- func_name = "abs"
-
- return var(func_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.0
-
- def map_node_coordinate_component(self, expr):
- return self.map_variable_ref_expr(
- expr, "grdg_ncc%d" % expr.axis)
-
- def map_common_subexpression(self, expr):
- raise ValueError("not expecting CSEs at this stage in the "
- "compilation process")
-
-
-# {{{ bessel handling
-
-class BesselFunction(ScalarCallable):
- def with_types(self, arg_id_to_dtype, clbl_inf_ctx):
- from loopy.types import NumpyType
-
- for i in arg_id_to_dtype:
- if not (-1 <= i <= 1):
- raise TypeError(f"{self.name} can only take 2 arguments.")
-
- if (arg_id_to_dtype.get(0) is None) or (arg_id_to_dtype.get(1) is None):
- # not enough info about input types
- return self, clbl_inf_ctx
-
- n_dtype = arg_id_to_dtype[0]
-
- # *technically* takes a float, but let's not worry about that.
- if n_dtype.numpy_dtype.kind != "i":
- raise TypeError(f"{self.name} expects an integer first argument")
-
- if self.name == "bessel_j":
- name_in_target = "bessel_jv"
- else:
- assert self.name == "bessel_y"
- name_in_target = "bessel_yn"
-
- return (self.copy(name_in_target=name_in_target,
- arg_id_to_dtype={-1: NumpyType(np.float64),
- 0: NumpyType(np.int32),
- 1: NumpyType(np.float64),
- }),
- clbl_inf_ctx)
-
- def generate_preambles(self, target):
- from loopy import PyOpenCLTarget
- if not isinstance(target, PyOpenCLTarget):
- raise NotImplementedError("Only the PyOpenCLTarget is supported as"
- "of now.")
-
- yield ("50-grudge-bessel", """//CL//
- #include <pyopencl-bessel-j.cl>
- #include <pyopencl-bessel-y.cl>
- """)
- return
-
-# }}}
-
-
-class ToLoopyInstructionMapper:
- def __init__(self, dd_inference_mapper):
- self.dd_inference_mapper = dd_inference_mapper
- self.function_registry = dd_inference_mapper.function_registry
- self.insn_count = 0
-
- def map_insn_assign(self, insn):
- from grudge.symbolic.primitives import OperatorBinding
-
- if (
- len(insn.exprs) == 1
- and (
- isinstance(insn.exprs[0], OperatorBinding)
- or is_external_call(
- insn.exprs[0], self.function_registry))):
- return insn
-
- # FIXME: These names and the size names could clash with user-given names.
- # Need better metadata tracking in loopy.
- iel = "iel"
- idof = "idof"
-
- temp_names = [
- name
- for name, dnr in zip(insn.names, insn.do_not_return)
- if dnr]
-
- from pymbolic import var
- expr_mapper = ToLoopyExpressionMapper(
- self.dd_inference_mapper, temp_names, (var(iel), var(idof)))
- 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.Optional(None) if dnr else lp.Optional(),
- no_sync_with=frozenset([
- ("*", "any"),
- ]),
- ))
-
- if not expr_mapper.non_scalar_vars:
- return insn
-
- knl = lp.make_kernel(
- "{[%(iel)s, %(idof)s]: "
- "0 <= %(iel)s < nelements and 0 <= %(idof)s < nunit_dofs}"
- % {"iel": iel, "idof": idof},
- insns,
-
- [
- lp.GlobalArg(name, shape=lp.auto, is_input=False)
- for name, dnr in zip(insn.names, insn.do_not_return)
- if not dnr
- ] + [...],
- name="grudge_assign_%d" % self.insn_count,
-
- # Single-insn kernels may have their no_sync_with resolve to an
- # empty set, that's OK.
- options=lp.Options(
- check_dep_resolution=False,
- return_dict=True,
- no_numpy=True,
- )
- )
-
- self.insn_count += 1
-
- from pytools import single_valued
- governing_dd = single_valued(
- self.dd_inference_mapper(expr)
- for expr in insn.exprs)
-
- knl = lp.register_callable(knl, "bessel_j", BesselFunction("bessel_j"))
- knl = lp.register_callable(knl, "bessel_y", BesselFunction("bessel_y"))
-
- input_mappings = {}
- output_mappings = {}
-
- from grudge.symbolic.mappers import DependencyMapper
- dep_mapper = DependencyMapper(composite_leaves=False)
-
- for expr, name in expr_mapper.expr_to_name.items():
- deps = dep_mapper(expr)
- assert len(deps) <= 1
- if not deps:
- is_output = False
- else:
- dep, = deps
- is_output = dep.name in insn.names
-
- if is_output:
- tgt_dict = output_mappings
- else:
- tgt_dict = input_mappings
-
- tgt_dict[name] = expr
-
- return LoopyKernelInstruction(
- LoopyKernelDescriptor(
- loopy_kernel=knl,
- input_mappings=input_mappings,
- output_mappings=output_mappings,
- fixed_arguments={},
- governing_dd=governing_dd)
- )
-
- def map_insn_rank_data_swap(self, insn):
- return insn
-
- 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):
- """
- .. attribute:: generating_discr_code
- """
-
- def get_expr_to_var(self, compiler):
- if self.generating_discr_code:
- return compiler.discr_expr_to_var
- else:
- return compiler.eval_expr_to_var
-
- def get_code_list(self, compiler):
- if self.generating_discr_code:
- return compiler.discr_code
- else:
- return compiler.eval_code
-
-
-class OperatorCompiler(mappers.IdentityMapper):
- def __init__(self, discr, function_registry,
- prefix="_expr", max_vectors_in_batch_expr=None):
- super().__init__()
- self.prefix = prefix
-
- self.max_vectors_in_batch_expr = max_vectors_in_batch_expr
-
- self.discr_code = []
- self.discr_scope_names_created = set()
- self.discr_scope_names_copied_to_eval = set()
- self.discr_expr_to_var = {}
-
- self.eval_code = []
- self.eval_expr_to_var = {}
-
- self.assigned_names = set()
-
- self.discr = discr
- self.function_registry = function_registry
-
- from pytools import UniqueNameGenerator
- self.name_gen = UniqueNameGenerator()
-
- # {{{ collect various optemplate components
-
- def collect_diff_ops(self, expr):
- return mappers.BoundOperatorCollector(sym.RefDiffOperatorBase)(expr)
-
- # }}}
-
- # {{{ top-level driver
-
- def __call__(self, expr):
- # Put the result expressions into variables as well.
- expr = sym.cse(expr, "_result")
-
- # from grudge.symbolic.mappers.type_inference import TypeInferrer
- # self.typedict = TypeInferrer()(expr)
-
- # Used for diff batching
- self.diff_ops = self.collect_diff_ops(expr)
-
- codegen_state = CodeGenerationState(generating_discr_code=False)
- # Finally, walk the expression and build the code.
- result = super().__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, self.function_registry)
-
- 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(discr_code,
- make_obj_array(
- [Variable(name)
- for name in self.discr_scope_names_copied_to_eval])),
- Code(eval_code, result))
-
- # }}}
-
- # {{{ variables and names
-
- def assign_to_new_var(self, codegen_state, expr, priority=0, prefix=None):
- # Observe that the only things that can be legally subscripted in
- # grudge are variables. All other expressions are broken down into
- # their scalar components.
- if isinstance(expr, (Variable, Subscript)):
- return expr
-
- new_name = self.name_gen(prefix if prefix is not None else "expr")
- codegen_state.get_code_list(self).append(Assign(
- (new_name,), (expr,), priority=priority))
-
- return Variable(new_name)
-
- # }}}
-
- # {{{ map_xxx routines
-
- def map_common_subexpression(self, expr, codegen_state):
- def get_rec_child(codegen_state):
- if isinstance(expr.child, sym.OperatorBinding):
- # We need to catch operator bindings here and
- # treat them specially. They get assigned to their
- # own variable by default, which would mean the
- # CSE prefix would be omitted, making the resulting
- # code less readable.
- return self.map_operator_binding(
- expr.child, codegen_state, name_hint=expr.prefix)
- else:
- return self.rec(expr.child, codegen_state)
-
- if expr.scope == sym.cse_scope.DISCRETIZATION:
- from pymbolic import var
- try:
- expr_name = self.discr._discr_scoped_subexpr_to_name[expr.child]
- except KeyError:
- expr_name = "discr." + self.discr._discr_scoped_name_gen(
- expr.prefix if expr.prefix is not None else "expr")
- self.discr._discr_scoped_subexpr_to_name[expr.child] = expr_name
-
- assert expr_name.startswith("discr.")
-
- priority = getattr(expr, "priority", 0)
-
- if expr_name not in self.discr_scope_names_created:
- new_codegen_state = codegen_state.copy(generating_discr_code=True)
- rec_child = get_rec_child(new_codegen_state)
-
- new_codegen_state.get_code_list(self).append(
- ToDiscretizationScopedAssign(
- (expr_name,), (rec_child,), priority=priority))
-
- self.discr_scope_names_created.add(expr_name)
-
- if codegen_state.generating_discr_code:
- return var(expr_name)
- else:
- if expr_name in self.discr_scope_names_copied_to_eval:
- return var(expr_name)
-
- self.eval_code.append(
- FromDiscretizationScopedAssign(
- expr_name, priority=priority))
-
- self.discr_scope_names_copied_to_eval.add(expr_name)
-
- return var(expr_name)
-
- else:
- expr_to_var = codegen_state.get_expr_to_var(self)
-
- try:
- return expr_to_var[expr.child]
- except KeyError:
- priority = getattr(expr, "priority", 0)
-
- rec_child = get_rec_child(codegen_state)
-
- cse_var = self.assign_to_new_var(
- codegen_state, rec_child,
- priority=priority, prefix=expr.prefix)
-
- expr_to_var[expr.child] = cse_var
- return cse_var
-
- def map_operator_binding(self, expr, codegen_state, name_hint=None):
- if isinstance(expr.op, sym.RefDiffOperatorBase):
- return self.map_ref_diff_op_binding(expr, codegen_state)
- elif isinstance(expr.op, sym.OppositePartitionFaceSwap):
- return self.map_rank_data_swap_binding(expr, codegen_state, name_hint)
- else:
- # make sure operator assignments stand alone and don't get muddled
- # up in vector math
- field_var = self.assign_to_new_var(
- codegen_state,
- self.rec(expr.field, codegen_state))
- result_var = self.assign_to_new_var(
- codegen_state,
- expr.op(field_var),
- prefix=name_hint)
- return result_var
-
- def map_call(self, expr, codegen_state):
- if is_function_loopyable(expr.function, self.function_registry):
- return super().map_call(expr, codegen_state)
- else:
- # If it's not a C-level function, it shouldn't get muddled up into
- # a vector math expression.
-
- return self.assign_to_new_var(
- codegen_state,
- type(expr)(
- expr.function,
- [self.assign_to_new_var(
- codegen_state,
- self.rec(par, codegen_state))
- for par in expr.parameters]))
-
- def map_ref_diff_op_binding(self, expr, codegen_state):
- expr_to_var = codegen_state.get_expr_to_var(self)
- try:
- return expr_to_var[expr]
- except KeyError:
- all_diffs = [diff
- for diff in self.diff_ops
- if diff.op.equal_except_for_axis(expr.op)
- and diff.field == expr.field]
-
- names = [self.name_gen("expr") for d in all_diffs]
-
- from pytools import single_valued
- op_class = single_valued(type(d.op) for d in all_diffs)
-
- codegen_state.get_code_list(self).append(
- DiffBatchAssign(
- names=names,
- op_class=op_class,
- operators=[d.op for d in all_diffs],
- field=self.rec(
- single_valued(d.field for d in all_diffs),
- codegen_state)))
-
- from pymbolic import var
- for n, d in zip(names, all_diffs):
- expr_to_var[d] = var(n)
-
- return expr_to_var[expr]
-
- def map_rank_data_swap_binding(self, expr, codegen_state, name_hint):
- expr_to_var = codegen_state.get_expr_to_var(self)
-
- try:
- return expr_to_var[expr]
- except KeyError:
- field = self.rec(expr.field, codegen_state)
- name = self.name_gen("raw_rank%02d_bdry_data" % expr.op.i_remote_part)
- field_insn = RankDataSwapAssign(name=name, field=field, op=expr.op)
- codegen_state.get_code_list(self).append(field_insn)
- field_var = Variable(field_insn.name)
- expr_to_var[expr] = self.assign_to_new_var(codegen_state,
- expr.op(field_var),
- prefix=name_hint)
- return expr_to_var[expr]
-
- # }}}
-
-# }}}
-
-
-# vim: foldmethod=marker
diff --git a/grudge/symbolic/dofdesc_inference.py b/grudge/symbolic/dofdesc_inference.py
deleted file mode 100644
index 6c0da95d..00000000
--- a/grudge/symbolic/dofdesc_inference.py
+++ /dev/null
@@ -1,232 +0,0 @@
-__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.dof_desc 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 = ""
-
- 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.discretization_tag != dd_b.discretization_tag:
- raise ValueError("mismatched discretization tags " + loc_str)
-
- return dd_a
-
-
-class InferrableMultiAssignment:
- """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, function_registry,
- 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 = {
- 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
-
- self.function_registry = function_registry
-
- 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
- map_max = map_sum
- map_min = 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.op.dd_in,
- str(expr)))
-
- return operator.dd_out
-
- def map_ones(self, expr):
- return expr.dd
-
- map_node_coordinate_component = map_ones
- map_signed_face_ones = map_ones
-
- def map_call(self, expr):
- arg_dds = [
- self.rec(par)
- for par in expr.parameters]
-
- return (
- self.function_registry[expr.function.name]
- .get_result_dofdesc(arg_dds))
-
- def map_common_subexpression_uncached(self, expr):
- return self.rec(expr.child)
-
- # }}}
-
- # {{{ instruction mappings
-
- def map_insn_assign(self, insn):
- return [
- (name, self.rec(expr))
- for name, expr in zip(insn.names, insn.exprs)
- ]
-
- def map_insn_rank_data_swap(self, insn):
- return [(insn.name, insn.dd_out)]
-
- 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
deleted file mode 100644
index 288cc30c..00000000
--- a/grudge/symbolic/mappers/__init__.py
+++ /dev/null
@@ -1,1298 +0,0 @@
-"""Mappers to transform symbolic operators."""
-
-__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 numpy as np
-import pymbolic.primitives
-import pymbolic.mapper.stringifier
-import pymbolic.mapper.evaluator
-import pymbolic.mapper.dependency
-import pymbolic.mapper.substitutor
-import pymbolic.mapper.constant_folder
-import pymbolic.mapper.constant_converter
-import pymbolic.mapper.flop_counter
-from pymbolic.mapper import CSECachingMapperMixin
-
-from grudge import sym
-import grudge.dof_desc as dof_desc
-import grudge.symbolic.operators as op
-from grudge.tools import OrderedSet
-
-
-# {{{ mixins
-
-class LocalOpReducerMixin:
- """Reduces calls to mapper methods for all local differentiation
- operators to a single mapper method, and likewise for mass
- operators.
- """
- # {{{ global differentiation
- def map_diff(self, expr, *args, **kwargs):
- return self.map_diff_base(expr, *args, **kwargs)
-
- def map_minv_st(self, expr, *args, **kwargs):
- return self.map_diff_base(expr, *args, **kwargs)
-
- def map_stiffness(self, expr, *args, **kwargs):
- return self.map_diff_base(expr, *args, **kwargs)
-
- def map_stiffness_t(self, expr, *args, **kwargs):
- return self.map_diff_base(expr, *args, **kwargs)
-
- def map_quad_stiffness_t(self, expr, *args, **kwargs):
- return self.map_diff_base(expr, *args, **kwargs)
- # }}}
-
- # {{{ global mass
- def map_mass_base(self, expr, *args, **kwargs):
- return self.map_elementwise_linear(expr, *args, **kwargs)
-
- def map_mass(self, expr, *args, **kwargs):
- return self.map_mass_base(expr, *args, **kwargs)
-
- def map_inverse_mass(self, expr, *args, **kwargs):
- return self.map_mass_base(expr, *args, **kwargs)
-
- def map_quad_mass(self, expr, *args, **kwargs):
- return self.map_mass_base(expr, *args, **kwargs)
- # }}}
-
- # {{{ reference differentiation
- def map_ref_diff(self, expr, *args, **kwargs):
- return self.map_ref_diff_base(expr, *args, **kwargs)
-
- def map_ref_stiffness_t(self, expr, *args, **kwargs):
- return self.map_ref_diff_base(expr, *args, **kwargs)
-
- def map_ref_quad_stiffness_t(self, expr, *args, **kwargs):
- return self.map_ref_diff_base(expr, *args, **kwargs)
- # }}}
-
- # {{{ reference mass
- def map_ref_mass_base(self, expr, *args, **kwargs):
- return self.map_elementwise_linear(expr, *args, **kwargs)
-
- def map_ref_mass(self, expr, *args, **kwargs):
- return self.map_ref_mass_base(expr, *args, **kwargs)
-
- def map_ref_inverse_mass(self, expr, *args, **kwargs):
- return self.map_ref_mass_base(expr, *args, **kwargs)
-
- def map_ref_quad_mass(self, expr, *args, **kwargs):
- return self.map_ref_mass_base(expr, *args, **kwargs)
- # }}}
-
-
-class FluxOpReducerMixin:
- """Reduces calls to mapper methods for all flux
- operators to a smaller number of mapper methods.
- """
- def map_flux(self, expr, *args, **kwargs):
- return self.map_flux_base(expr, *args, **kwargs)
-
- def map_bdry_flux(self, expr, *args, **kwargs):
- return self.map_flux_base(expr, *args, **kwargs)
-
- def map_quad_flux(self, expr, *args, **kwargs):
- return self.map_flux_base(expr, *args, **kwargs)
-
- def map_quad_bdry_flux(self, expr, *args, **kwargs):
- return self.map_flux_base(expr, *args, **kwargs)
-
-
-class OperatorReducerMixin(LocalOpReducerMixin, FluxOpReducerMixin):
- """Reduces calls to *any* operator mapping function to just one."""
- def _map_op_base(self, expr, *args, **kwargs):
- return self.map_operator(expr, *args, **kwargs)
-
- map_elementwise_linear = _map_op_base
-
- map_projection = _map_op_base
-
- map_elementwise_sum = _map_op_base
- map_elementwise_min = _map_op_base
- map_elementwise_max = _map_op_base
-
- map_nodal_sum = _map_op_base
- map_nodal_min = _map_op_base
- map_nodal_max = _map_op_base
-
- map_stiffness = _map_op_base
- map_diff = _map_op_base
- map_stiffness_t = _map_op_base
-
- map_ref_diff = _map_op_base
- map_ref_stiffness_t = _map_op_base
-
- map_mass = _map_op_base
- map_inverse_mass = _map_op_base
- map_ref_mass = _map_op_base
- map_ref_inverse_mass = _map_op_base
-
- map_opposite_partition_face_swap = _map_op_base
- map_opposite_interior_face_swap = _map_op_base
- map_face_mass_operator = _map_op_base
- map_ref_face_mass_operator = _map_op_base
-
-
-class CombineMapperMixin:
- def map_operator_binding(self, expr):
- return self.combine([self.rec(expr.op), self.rec(expr.field)])
-
-
-class IdentityMapperMixin(LocalOpReducerMixin, FluxOpReducerMixin):
- def map_operator_binding(self, expr, *args, **kwargs):
- assert not isinstance(self, BoundOpMapperMixin), \
- "IdentityMapper instances cannot be combined with " \
- "the BoundOpMapperMixin"
-
- return type(expr)(
- self.rec(expr.op, *args, **kwargs),
- self.rec(expr.field, *args, **kwargs))
-
- # {{{ operators
-
- def map_elementwise_linear(self, expr, *args, **kwargs):
- assert not isinstance(self, BoundOpMapperMixin), \
- "IdentityMapper instances cannot be combined with " \
- "the BoundOpMapperMixin"
-
- # it's a leaf--no changing children
- return expr
-
- map_projection = map_elementwise_linear
-
- map_elementwise_sum = map_elementwise_linear
- map_elementwise_min = map_elementwise_linear
- map_elementwise_max = map_elementwise_linear
-
- map_nodal_sum = map_elementwise_linear
- map_nodal_min = map_elementwise_linear
- map_nodal_max = map_elementwise_linear
-
- map_stiffness = map_elementwise_linear
- map_diff = map_elementwise_linear
- map_stiffness_t = map_elementwise_linear
-
- map_ref_diff = map_elementwise_linear
- map_ref_stiffness_t = map_elementwise_linear
-
- map_mass = map_elementwise_linear
- map_inverse_mass = map_elementwise_linear
- map_ref_mass = map_elementwise_linear
- map_ref_inverse_mass = map_elementwise_linear
-
- map_opposite_partition_face_swap = map_elementwise_linear
- map_opposite_interior_face_swap = map_elementwise_linear
- map_face_mass_operator = map_elementwise_linear
- map_ref_face_mass_operator = map_elementwise_linear
-
- # }}}
-
- # {{{ primitives
-
- def map_grudge_variable(self, expr, *args, **kwargs):
- # it's a leaf--no changing children
- return expr
-
- map_function_symbol = map_grudge_variable
- map_ones = map_grudge_variable
- map_signed_face_ones = map_grudge_variable
- map_node_coordinate_component = map_grudge_variable
-
- # }}}
-
-
-class BoundOpMapperMixin:
- def map_operator_binding(self, expr, *args, **kwargs):
- return getattr(self, expr.op.mapper_method)(
- expr.op, expr.field, *args, **kwargs)
-
-# }}}
-
-
-# {{{ basic mappers
-
-class CombineMapper(CombineMapperMixin, pymbolic.mapper.CombineMapper):
- pass
-
-
-class DependencyMapper(
- CombineMapperMixin,
- pymbolic.mapper.dependency.DependencyMapper,
- OperatorReducerMixin):
- def __init__(self,
- include_operator_bindings=True,
- composite_leaves=None,
- **kwargs):
- if composite_leaves is False:
- include_operator_bindings = False
- if composite_leaves is True:
- include_operator_bindings = True
-
- pymbolic.mapper.dependency.DependencyMapper.__init__(self,
- composite_leaves=composite_leaves, **kwargs)
-
- self.include_operator_bindings = include_operator_bindings
-
- def map_operator_binding(self, expr):
- if self.include_operator_bindings:
- return {expr}
- else:
- return CombineMapperMixin.map_operator_binding(self, expr)
-
- def map_operator(self, expr):
- return set()
-
- def map_grudge_variable(self, expr):
- return {expr}
-
- def _map_leaf(self, expr):
- return set()
-
- map_ones = _map_leaf
- map_signed_face_ones = _map_leaf
- map_node_coordinate_component = _map_leaf
-
-
-class FlopCounter(
- CombineMapperMixin,
- pymbolic.mapper.flop_counter.FlopCounter):
- def map_operator_binding(self, expr):
- return self.rec(expr.field)
-
- def map_grudge_variable(self, expr):
- return 0
-
- def map_function_symbol(self, expr):
- return 1
-
- def map_ones(self, expr):
- return 0
-
- def map_node_coordinate_component(self, expr):
- return 0
-
-
-class IdentityMapper(
- IdentityMapperMixin,
- pymbolic.mapper.IdentityMapper):
- pass
-
-
-class SubstitutionMapper(pymbolic.mapper.substitutor.SubstitutionMapper,
- IdentityMapperMixin):
- pass
-
-
-class CSERemover(IdentityMapper):
- def map_common_subexpression(self, expr):
- return self.rec(expr.child)
-
-# }}}
-
-
-# {{{ operator binder
-
-class OperatorBinder(CSECachingMapperMixin, IdentityMapper):
- map_common_subexpression_uncached = \
- IdentityMapper.map_common_subexpression
-
- def map_product(self, expr):
- if len(expr.children) == 0:
- return expr
-
- from pymbolic.primitives import flattened_product, Product
-
- first = expr.children[0]
- if isinstance(first, op.Operator):
- prod = flattened_product(expr.children[1:])
- if isinstance(prod, Product) and len(prod.children) > 1:
- from warnings import warn
- warn("Binding '%s' to more than one "
- "operand in a product is ambiguous - "
- "use the parenthesized form instead."
- % first)
- return sym.OperatorBinding(first, self.rec(prod))
- else:
- return self.rec(first) * self.rec(flattened_product(expr.children[1:]))
-
-# }}}
-
-
-# {{{ dof desc (dd) replacement
-
-class DOFDescReplacer(IdentityMapper):
- def __init__(self, prev_dd, new_dd):
- self.prev_dd = prev_dd
- self.new_dd = new_dd
-
- def map_operator_binding(self, expr):
- if (isinstance(expr.op, op.OppositeInteriorFaceSwap)
- and expr.op.dd_in == self.prev_dd
- and expr.op.dd_out == self.prev_dd):
- field = self.rec(expr.field)
- return op.OppositePartitionFaceSwap(dd_in=self.new_dd,
- dd_out=self.new_dd)(field)
- elif (isinstance(expr.op, op.ProjectionOperator)
- and expr.op.dd_out == self.prev_dd):
- return op.ProjectionOperator(dd_in=expr.op.dd_in,
- dd_out=self.new_dd)(expr.field)
- elif (isinstance(expr.op, op.RefDiffOperatorBase)
- and expr.op.dd_out == self.prev_dd
- and expr.op.dd_in == self.prev_dd):
- return type(expr.op)(expr.op.rst_axis,
- dd_in=self.new_dd,
- dd_out=self.new_dd)(self.rec(expr.field))
-
- def map_node_coordinate_component(self, expr):
- if expr.dd == self.prev_dd:
- return type(expr)(expr.axis, self.new_dd)
-
-# }}}
-
-
-# {{{ mappers for distributed computation
-
-class OppositeInteriorFaceSwapUniqueIDAssigner(
- CSECachingMapperMixin, IdentityMapper):
- map_common_subexpression_uncached = IdentityMapper.map_common_subexpression
-
- def __init__(self):
- super().__init__()
- self._next_id = 0
- self.seen_ids = set()
-
- def next_id(self):
- while self._next_id in self.seen_ids:
- self._next_id += 1
-
- result = self._next_id
- self._next_id += 1
- self.seen_ids.add(result)
-
- return result
-
- def map_opposite_interior_face_swap(self, expr):
- if expr.unique_id is not None:
- if expr.unique_id in self.seen_ids:
- raise ValueError("OppositeInteriorFaceSwap unique ID '%d' "
- "is not unique" % expr.unique_id)
-
- self.seen_ids.add(expr.unique_id)
- return expr
-
- else:
- return type(expr)(expr.dd_in, expr.dd_out, self.next_id())
-
-
-class DistributedMapper(CSECachingMapperMixin, IdentityMapper):
- map_common_subexpression_uncached = IdentityMapper.map_common_subexpression
-
- def __init__(self, connected_parts):
- self.connected_parts = connected_parts
-
- def map_operator_binding(self, expr):
- from meshmode.mesh import BTAG_PARTITION
- from meshmode.discretization.connection import (FACE_RESTR_ALL,
- FACE_RESTR_INTERIOR)
- if (isinstance(expr.op, op.ProjectionOperator)
- and expr.op.dd_in.domain_tag is FACE_RESTR_INTERIOR
- and expr.op.dd_out.domain_tag is FACE_RESTR_ALL):
- distributed_work = 0
- for i_remote_part in self.connected_parts:
- mapped_field = RankGeometryChanger(i_remote_part)(expr.field)
- btag_part = BTAG_PARTITION(i_remote_part)
- distributed_work += op.ProjectionOperator(dd_in=btag_part,
- dd_out=expr.op.dd_out)(mapped_field)
- return expr + distributed_work
- else:
- return IdentityMapper.map_operator_binding(self, expr)
-
-
-class RankGeometryChanger(CSECachingMapperMixin, IdentityMapper):
- map_common_subexpression_uncached = IdentityMapper.map_common_subexpression
-
- def __init__(self, i_remote_part):
- from meshmode.discretization.connection import FACE_RESTR_INTERIOR
- from meshmode.mesh import BTAG_PARTITION
-
- self.prev_dd = dof_desc.as_dofdesc(FACE_RESTR_INTERIOR)
- self.new_dd = dof_desc.as_dofdesc(BTAG_PARTITION(i_remote_part))
-
- def _raise_unable(self, expr):
- raise ValueError("encountered '%s' in updating subexpression for "
- "changed geometry (likely for distributed computation); "
- "unable to adapt from '%s' to '%s'"
- % (str(expr), self.prev_dd, self.new_dd))
-
- def map_operator_binding(self, expr):
- if (isinstance(expr.op, op.OppositeInteriorFaceSwap)
- and expr.op.dd_in == self.prev_dd
- and expr.op.dd_out == self.prev_dd):
- field = self.rec(expr.field)
- return op.OppositePartitionFaceSwap(
- dd_in=self.new_dd,
- dd_out=self.new_dd,
- unique_id=expr.op.unique_id)(field)
- elif (isinstance(expr.op, op.ProjectionOperator)
- and expr.op.dd_out == self.prev_dd):
- return op.ProjectionOperator(dd_in=expr.op.dd_in,
- dd_out=self.new_dd)(expr.field)
- elif (isinstance(expr.op, op.RefDiffOperator)
- and expr.op.dd_out == self.prev_dd
- and expr.op.dd_in == self.prev_dd):
- return op.RefDiffOperator(expr.op.rst_axis,
- dd_in=self.new_dd,
- dd_out=self.new_dd)(self.rec(expr.field))
- else:
- self._raise_unable(expr)
-
- def map_grudge_variable(self, expr):
- self._raise_unable(expr)
-
- def map_node_coordinate_component(self, expr):
- if expr.dd == self.prev_dd:
- return type(expr)(expr.axis, self.new_dd)
- else:
- self._raise_unable(expr)
-
-# }}}
-
-
-# {{{ operator specializer
-
-class OperatorSpecializer(CSECachingMapperMixin, IdentityMapper):
- """Guided by a typedict obtained through type inference (i.e. by
- :class:`grudge.symbolic.mappers.type_inference.TypeInferrrer`),
- substitutes more specialized operators for generic ones.
-
- For example, if type inference has determined that a differentiation
- operator is applied to an argument on a quadrature grid, this
- differentiation operator is then swapped out for a *quadrature*
- differentiation operator.
- """
-
- def __init__(self, typedict):
- """
- :param typedict: generated by
- :class:`grudge.symbolic.mappers.type_inference.TypeInferrer`.
- """
- self.typedict = typedict
-
- map_common_subexpression_uncached = \
- IdentityMapper.map_common_subexpression
-
- def map_operator_binding(self, expr):
- from grudge.symbolic.primitives import TracePair
-
- from grudge.symbolic.mappers.type_inference import (
- type_info, QuadratureRepresentation)
-
- # {{{ figure out field type
- try:
- field_type = self.typedict[expr.field]
- except TypeError:
- # numpy arrays are not hashable
- # has_quad_operand remains unset
-
- assert isinstance(expr.field, np.ndarray)
- else:
- try:
- field_repr_tag = field_type.repr_tag
- except AttributeError:
- # boundary pairs are not assigned types
- assert isinstance(expr.field, TracePair)
- has_quad_operand = False
- else:
- has_quad_operand = isinstance(field_repr_tag,
- QuadratureRepresentation)
- # }}}
-
- # 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, op.MassOperator) and has_quad_operand:
- return op.QuadratureMassOperator(
- field_repr_tag.discretization_tag)(self.rec(expr.field))
-
- elif isinstance(expr.op, op.RefMassOperator) and has_quad_operand:
- return op.RefQuadratureMassOperator(
- field_repr_tag.discretization_tag)(self.rec(expr.field))
-
- elif (isinstance(expr.op, op.StiffnessTOperator) and has_quad_operand):
- return op.QuadratureStiffnessTOperator(
- expr.op.xyz_axis, field_repr_tag.discretization_tag)(
- self.rec(expr.field))
-
- elif (isinstance(expr.op, op.RefStiffnessTOperator)
- and has_quad_operand):
- return op.RefQuadratureStiffnessTOperator(
- expr.op.xyz_axis, field_repr_tag.discretization_tag)(
- self.rec(expr.field))
-
- elif (isinstance(expr.op, op.QuadratureGridUpsampler)
- and isinstance(field_type, type_info.BoundaryVectorBase)):
- # potential shortcut:
- # if (isinstance(expr.field, OperatorBinding)
- # and isinstance(expr.field.op, RestrictToBoundary)):
- # return QuadratureRestrictToBoundary(
- # expr.field.op.tag, expr.op.discretization_tag)(
- # self.rec(expr.field.field))
-
- return op.QuadratureBoundaryGridUpsampler(
- expr.op.discretization_tag, field_type.boundary_tag)(expr.field)
- # }}}
-
- elif isinstance(expr.op, op.RestrictToBoundary) and has_quad_operand:
- raise TypeError("RestrictToBoundary cannot be applied to "
- "quadrature-based operands--use QuadUpsample(Boundarize(...))")
-
- else:
- return IdentityMapper.map_operator_binding(self, expr)
-
-# }}}
-
-
-# {{{ global-to-reference mapper
-
-class GlobalToReferenceMapper(CSECachingMapperMixin, IdentityMapper):
- """Maps operators that apply on the global function space down to operators on
- reference elements, together with explicit multiplication by geometric factors.
- """
-
- def __init__(self, dcoll):
- CSECachingMapperMixin.__init__(self)
- IdentityMapper.__init__(self)
-
- self.ambient_dim = dcoll.ambient_dim
- self.dim = dcoll.dim
-
- map_common_subexpression_uncached = \
- IdentityMapper.map_common_subexpression
-
- def map_operator_binding(self, expr):
- # Global-to-reference is run after operator specialization, so
- # if we encounter non-quadrature operators here, we know they
- # must be nodal.
-
- dd_in = expr.op.dd_in
- dd_out = expr.op.dd_out
-
- if dd_in.is_volume():
- dim = self.dim
- else:
- dim = self.dim - 1
-
- jac_in = sym.area_element(self.ambient_dim, dim, dd=dd_in)
- jac_noquad = sym.area_element(self.ambient_dim, dim,
- dd=dd_in.with_discr_tag(dof_desc.DISCR_TAG_BASE))
-
- def rewrite_derivative(ref_class, field, dd_in, with_jacobian=True):
- def imd(rst):
- return sym.inverse_surface_metric_derivative(
- rst, expr.op.xyz_axis,
- ambient_dim=self.ambient_dim, dim=self.dim,
- dd=dd_in)
-
- rec_field = self.rec(field)
- if with_jacobian:
- jac_tag = sym.area_element(self.ambient_dim, self.dim, dd=dd_in)
- rec_field = jac_tag * rec_field
-
- return sum(
- ref_class(rst_axis, dd_in=dd_in)(rec_field * imd(rst_axis))
- for rst_axis in range(self.dim))
- else:
- return sum(
- ref_class(rst_axis, dd_in=dd_in)(rec_field) * imd(rst_axis)
- for rst_axis in range(self.dim))
-
- if isinstance(expr.op, op.MassOperator):
- return op.RefMassOperator(dd_in, dd_out)(
- jac_in * self.rec(expr.field))
-
- elif isinstance(expr.op, op.InverseMassOperator):
- # based on https://arxiv.org/pdf/1608.03836.pdf
- return (
- 1.0/jac_in * op.RefInverseMassOperator(dd_in, dd_out)(
- self.rec(expr.field)))
-
- elif isinstance(expr.op, op.FaceMassOperator):
- jac_in_surf = sym.area_element(
- self.ambient_dim, self.dim - 1, dd=dd_in)
- return op.RefFaceMassOperator(dd_in, dd_out)(
- jac_in_surf * self.rec(expr.field))
-
- elif isinstance(expr.op, op.StiffnessOperator):
- return op.RefMassOperator(dd_in=dd_in, dd_out=dd_out)(
- jac_noquad
- * self.rec(
- op.DiffOperator(expr.op.xyz_axis)(expr.field)))
-
- elif isinstance(expr.op, op.DiffOperator):
- return rewrite_derivative(
- op.RefDiffOperator,
- expr.field, dd_in=dd_in, with_jacobian=False)
-
- elif isinstance(expr.op, op.StiffnessTOperator):
- return rewrite_derivative(
- op.RefStiffnessTOperator,
- expr.field, dd_in=dd_in)
-
- elif isinstance(expr.op, op.MInvSTOperator):
- return self.rec(
- op.InverseMassOperator()(
- op.StiffnessTOperator(expr.op.xyz_axis)(
- self.rec(expr.field))))
-
- else:
- return IdentityMapper.map_operator_binding(self, expr)
-
-# }}}
-
-
-# {{{ stringification ---------------------------------------------------------
-
-class StringifyMapper(pymbolic.mapper.stringifier.StringifyMapper):
- def _format_dd(self, dd):
- def fmt(s):
- if isinstance(s, type):
- return s.__name__
- else:
- return repr(s)
-
- from meshmode.mesh import BTAG_PARTITION
- from meshmode.discretization.connection import (
- FACE_RESTR_ALL, FACE_RESTR_INTERIOR)
- if dd.domain_tag is None:
- result = "?"
- elif dd.domain_tag is dof_desc.DTAG_VOLUME_ALL:
- result = "vol"
- elif dd.domain_tag is dof_desc.DTAG_SCALAR:
- result = "scalar"
- elif dd.domain_tag is FACE_RESTR_ALL:
- result = "all_faces"
- elif dd.domain_tag is FACE_RESTR_INTERIOR:
- result = "int_faces"
- elif isinstance(dd.domain_tag, BTAG_PARTITION):
- result = "part%d_faces" % dd.domain_tag.part_nr
- else:
- result = fmt(dd.domain_tag)
-
- if dd.discretization_tag is None:
- pass
- elif dd.discretization_tag is dof_desc.DISCR_TAG_BASE:
- result += "q"
- else:
- result += "Q"+fmt(dd.discretization_tag)
-
- return result
-
- def _format_op_dd(self, op):
- return ":{}->{}".format(
- self._format_dd(op.dd_in),
- self._format_dd(op.dd_out))
-
- # {{{ elementwise ops
-
- def map_elementwise_sum(self, expr, enclosing_prec):
- return "ElementwiseSum" + self._format_op_dd(expr)
-
- def map_elementwise_max(self, expr, enclosing_prec):
- return "ElementwiseMax" + self._format_op_dd(expr)
-
- def map_elementwise_min(self, expr, enclosing_prec):
- return "ElementwiseMin" + self._format_op_dd(expr)
-
- # }}}
-
- # {{{ nodal ops
-
- def map_nodal_sum(self, expr, enclosing_prec):
- return "NodalSum" + self._format_op_dd(expr)
-
- def map_nodal_max(self, expr, enclosing_prec):
- return "NodalMax" + self._format_op_dd(expr)
-
- def map_nodal_min(self, expr, enclosing_prec):
- return "NodalMin" + self._format_op_dd(expr)
-
- # }}}
-
- # {{{ global differentiation
-
- def map_diff(self, expr, enclosing_prec):
- return "Diffx%d%s" % (expr.xyz_axis, self._format_op_dd(expr))
-
- def map_minv_st(self, expr, enclosing_prec):
- return "MInvSTx%d%s" % (expr.xyz_axis, self._format_op_dd(expr))
-
- def map_stiffness(self, expr, enclosing_prec):
- return "Stiffx%d%s" % (expr.xyz_axis, self._format_op_dd(expr))
-
- def map_stiffness_t(self, expr, enclosing_prec):
- return "StiffTx%d%s" % (expr.xyz_axis, self._format_op_dd(expr))
-
- # }}}
-
- # {{{ global mass
-
- def map_mass(self, expr, enclosing_prec):
- return "M"
-
- def map_inverse_mass(self, expr, enclosing_prec):
- return "InvM"
-
- # }}}
-
- # {{{ reference differentiation
-
- def map_ref_diff(self, expr, enclosing_prec):
- return "Diffr%d%s" % (expr.rst_axis, self._format_op_dd(expr))
-
- def map_ref_stiffness_t(self, expr, enclosing_prec):
- return "StiffTr%d%s" % (expr.rst_axis, self._format_op_dd(expr))
-
- # }}}
-
- # {{{ reference mass
-
- def map_ref_mass(self, expr, enclosing_prec):
- return "RefM" + self._format_op_dd(expr)
-
- def map_ref_inverse_mass(self, expr, enclosing_prec):
- return "RefInvM" + self._format_op_dd(expr)
-
- # }}}
-
- def map_elementwise_linear(self, expr, enclosing_prec):
- return "ElWLin:{}{}".format(
- expr.__class__.__name__,
- self._format_op_dd(expr))
-
- # {{{ flux
-
- def map_face_mass_operator(self, expr, enclosing_prec):
- return "FaceM" + self._format_op_dd(expr)
-
- def map_ref_face_mass_operator(self, expr, enclosing_prec):
- return "RefFaceM" + self._format_op_dd(expr)
-
- def map_opposite_partition_face_swap(self, expr, enclosing_prec):
- return "PartSwap" + self._format_op_dd(expr)
-
- def map_opposite_interior_face_swap(self, expr, enclosing_prec):
- return "OppSwap" + self._format_op_dd(expr)
-
- # }}}
-
- def map_ones(self, expr, enclosing_prec):
- return "Ones:" + self._format_dd(expr.dd)
-
- def map_signed_face_ones(self, expr, enclosing_prec):
- return "SignedOnes:" + self._format_dd(expr.dd)
-
- # {{{ geometry data
-
- def map_node_coordinate_component(self, expr, enclosing_prec):
- return "x[%d]@%s" % (expr.axis, self._format_dd(expr.dd))
-
- # }}}
-
- def map_operator_binding(self, expr, enclosing_prec):
- from pymbolic.mapper.stringifier import PREC_NONE
- return "<{}>({})".format(
- self.rec(expr.op, PREC_NONE),
- self.rec(expr.field, PREC_NONE))
-
- def map_grudge_variable(self, expr, enclosing_prec):
- return "{}:{}".format(expr.name, self._format_dd(expr.dd))
-
- def map_function_symbol(self, expr, enclosing_prec):
- return expr.name
-
- def map_projection(self, expr, enclosing_prec):
- return "Project" + self._format_op_dd(expr)
-
-
-class PrettyStringifyMapper(
- pymbolic.mapper.stringifier.CSESplittingStringifyMapperMixin,
- StringifyMapper):
- pass
-
-
-class NoCSEStringifyMapper(StringifyMapper):
- def map_common_subexpression(self, expr, enclosing_prec):
- return self.rec(expr.child, enclosing_prec)
-
-# }}}
-
-
-# {{{ quadrature support
-
-class QuadratureCheckerAndRemover(CSECachingMapperMixin, IdentityMapper):
- """Checks whether all quadratu
- """
-
- def __init__(self, discr_tag_to_group_factory):
- IdentityMapper.__init__(self)
- CSECachingMapperMixin.__init__(self)
- self.discr_tag_to_group_factory = discr_tag_to_group_factory
-
- map_common_subexpression_uncached = \
- IdentityMapper.map_common_subexpression
-
- def _process_dd(self, dd, location_descr):
-
- if dd.discretization_tag is not dof_desc.DISCR_TAG_BASE:
- if dd.discretization_tag not in self.discr_tag_to_group_factory:
- raise ValueError("found unknown quadrature tag '%s' in '%s'"
- % (dd.discretization_tag, location_descr))
-
- grp_factory = self.discr_tag_to_group_factory[dd.discretization_tag]
- if grp_factory is None:
- dd = dof_desc.DOFDesc(dd.domain_tag, dof_desc.DISCR_TAG_BASE)
-
- return dd
-
- def map_operator_binding(self, expr):
- dd_in_orig = dd_in = expr.op.dd_in
- dd_out_orig = dd_out = expr.op.dd_out
-
- dd_in = self._process_dd(dd_in, "dd_in of %s" % type(expr.op).__name__)
- dd_out = self._process_dd(dd_out, "dd_out of %s" % type(expr.op).__name__)
-
- if dd_in_orig == dd_in and dd_out_orig == dd_out:
- # unchanged
- return IdentityMapper.map_operator_binding(self, expr)
-
- import grudge.symbolic.operators as op
- # changed
-
- if dd_in == dd_out and isinstance(expr.op, op.ProjectionOperator):
- # This used to be to-quad interpolation and has become a no-op.
- # Remove it.
- return self.rec(expr.field)
-
- if isinstance(expr.op, op.StiffnessTOperator):
- new_op = type(expr.op)(expr.op.xyz_axis, dd_in, dd_out)
- elif isinstance(expr.op, (op.FaceMassOperator, op.ProjectionOperator)):
- new_op = type(expr.op)(dd_in, dd_out)
- else:
- raise NotImplementedError("do not know how to modify dd_in and dd_out "
- "in %s" % type(expr.op).__name__)
-
- return new_op(self.rec(expr.field))
-
- def map_ones(self, expr):
- dd = self._process_dd(expr.dd, location_descr=type(expr).__name__)
- return type(expr)(dd)
-
- def map_node_coordinate_component(self, expr):
- dd = self._process_dd(expr.dd, location_descr=type(expr).__name__)
- return type(expr)(expr.axis, dd)
-
-# }}}
-
-
-# {{{ simplification / optimization
-
-class ConstantToNumpyConversionMapper(
- CSECachingMapperMixin,
- pymbolic.mapper.constant_converter.ConstantToNumpyConversionMapper,
- IdentityMapperMixin):
- map_common_subexpression_uncached = (
- pymbolic.mapper.constant_converter
- .ConstantToNumpyConversionMapper
- .map_common_subexpression)
-
-
-class CommutativeConstantFoldingMapper(
- pymbolic.mapper.constant_folder.CommutativeConstantFoldingMapper,
- IdentityMapperMixin):
-
- def __init__(self):
- pymbolic.mapper.constant_folder\
- .CommutativeConstantFoldingMapper.__init__(self)
- self.dep_mapper = DependencyMapper()
-
- def is_constant(self, expr):
- return not bool(self.dep_mapper(expr))
-
- def map_operator_binding(self, expr):
- field = self.rec(expr.field)
-
- from grudge.tools import is_zero
- if is_zero(field):
- return 0
-
- return expr.op(field)
-
-
-class EmptyFluxKiller(CSECachingMapperMixin, IdentityMapper):
- def __init__(self, mesh):
- IdentityMapper.__init__(self)
- self.mesh = mesh
-
- map_common_subexpression_uncached = \
- IdentityMapper.map_common_subexpression
-
- def map_operator_binding(self, expr):
- from meshmode.mesh import is_boundary_tag_empty
- if (isinstance(expr.op, sym.ProjectionOperator)
- and expr.op.dd_out.is_boundary_or_partition_interface()):
- domain_tag = expr.op.dd_out.domain_tag
- assert isinstance(domain_tag, dof_desc.DTAG_BOUNDARY)
- if is_boundary_tag_empty(self.mesh, domain_tag.tag):
- return 0
-
- return IdentityMapper.map_operator_binding(self, expr)
-
-
-class _InnerDerivativeJoiner(pymbolic.mapper.RecursiveMapper):
- def map_operator_binding(self, expr, derivatives):
- if isinstance(expr.op, op.DifferentiationOperator):
- derivatives.setdefault(expr.op, []).append(expr.field)
- return 0
- else:
- return DerivativeJoiner()(expr)
-
- def map_common_subexpression(self, expr, derivatives):
- # no use preserving these if we're moving derivatives around
- return self.rec(expr.child, derivatives)
-
- def map_sum(self, expr, derivatives):
- from pymbolic.primitives import flattened_sum
- return flattened_sum(tuple(
- self.rec(child, derivatives) for child in expr.children))
-
- def map_product(self, expr, derivatives):
- from grudge.symbolic.tools import is_scalar
- from pytools import partition
- scalars, nonscalars = partition(is_scalar, expr.children)
-
- if len(nonscalars) != 1:
- return DerivativeJoiner()(expr)
- else:
- from pymbolic import flattened_product
- factor = flattened_product(scalars)
- nonscalar, = nonscalars
-
- sub_derivatives = {}
- nonscalar = self.rec(nonscalar, sub_derivatives)
-
- def do_map(expr):
- if is_scalar(expr):
- return expr
- else:
- return self.rec(expr, derivatives)
-
- for operator, operands in sub_derivatives.items():
- for operand in operands:
- derivatives.setdefault(operator, []).append(
- factor*operand)
-
- return factor*nonscalar
-
- def map_constant(self, expr, *args):
- return DerivativeJoiner()(expr)
-
- def map_scalar_parameter(self, expr, *args):
- return DerivativeJoiner()(expr)
-
- def map_if(self, expr, *args):
- return DerivativeJoiner()(expr)
-
- def map_power(self, expr, *args):
- return DerivativeJoiner()(expr)
-
- # these two are necessary because they're forwarding targets
- def map_algebraic_leaf(self, expr, *args):
- return DerivativeJoiner()(expr)
-
- def map_quotient(self, expr, *args):
- return DerivativeJoiner()(expr)
-
- map_node_coordinate_component = map_algebraic_leaf
-
-
-class DerivativeJoiner(CSECachingMapperMixin, IdentityMapper):
- r"""Joins derivatives:
-
- .. math::
-
- \frac{\partial A}{\partial x} + \frac{\partial B}{\partial x}
- \rightarrow
- \frac{\partial (A+B)}{\partial x}.
- """
- map_common_subexpression_uncached = \
- IdentityMapper.map_common_subexpression
-
- def map_sum(self, expr):
- idj = _InnerDerivativeJoiner()
-
- def invoke_idj(expr):
- sub_derivatives = {}
- result = idj(expr, sub_derivatives)
- if not sub_derivatives:
- return expr
- else:
- for operator, operands in sub_derivatives.items():
- derivatives.setdefault(operator, []).extend(operands)
-
- return result
-
- derivatives = {}
- new_children = [invoke_idj(child)
- for child in expr.children]
-
- for operator, operands in derivatives.items():
- new_children.insert(0, operator(
- sum(self.rec(operand) for operand in operands)))
-
- from pymbolic.primitives import flattened_sum
- return flattened_sum(new_children)
-
-
-class _InnerInverseMassContractor(pymbolic.mapper.RecursiveMapper):
- def __init__(self, outer_mass_contractor):
- self.outer_mass_contractor = outer_mass_contractor
- self.extra_operator_count = 0
-
- def map_constant(self, expr):
- from grudge.tools import is_zero
-
- if is_zero(expr):
- return 0
- else:
- return op.OperatorBinding(
- op.InverseMassOperator(),
- self.outer_mass_contractor(expr))
-
- def map_algebraic_leaf(self, expr):
- return op.OperatorBinding(
- op.InverseMassOperator(),
- self.outer_mass_contractor(expr))
-
- def map_operator_binding(self, binding):
- if isinstance(binding.op, op.MassOperator):
- return binding.field
- elif isinstance(binding.op, op.StiffnessOperator):
- return op.DifferentiationOperator(binding.op.xyz_axis)(
- self.outer_mass_contractor(binding.field))
- elif isinstance(binding.op, op.StiffnessTOperator):
- return op.MInvSTOperator(binding.op.xyz_axis)(
- self.outer_mass_contractor(binding.field))
- elif isinstance(binding.op, op.FluxOperator):
- assert not binding.op.is_lift
-
- return op.FluxOperator(binding.op.flux, is_lift=True)(
- self.outer_mass_contractor(binding.field))
- elif isinstance(binding.op, op.BoundaryFluxOperator):
- assert not binding.op.is_lift
-
- 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 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):
- def is_scalar(expr):
- return isinstance(expr, (int, float, complex, sym.ScalarParameter))
-
- from pytools import len_iterable
- nonscalar_count = len_iterable(ch
- for ch in expr.children
- if not is_scalar(ch))
-
- if nonscalar_count > 1:
- # too complicated, don't touch it
- self.extra_operator_count += 1
- return op.InverseMassOperator()(
- self.outer_mass_contractor(expr))
- else:
- def do_map(expr):
- if is_scalar(expr):
- return expr
- else:
- return self.rec(expr)
- return expr.__class__(tuple(
- do_map(child) for child in expr.children))
-
-
-class InverseMassContractor(CSECachingMapperMixin, IdentityMapper):
- # assumes all operators to be bound
- map_common_subexpression_uncached = \
- IdentityMapper.map_common_subexpression
-
- def map_operator_binding(self, binding):
- # we only care about bindings of inverse mass operators
-
- if isinstance(binding.op, op.InverseMassOperator):
- iimc = _InnerInverseMassContractor(self)
- proposed_result = iimc(binding.field)
- if iimc.extra_operator_count > 1:
- # We're introducing more work than we're saving.
- # Don't perform the simplification
- return binding.op(self.rec(binding.field))
- else:
- return proposed_result
- else:
- return binding.op(self.rec(binding.field))
-
-# }}}
-
-
-# {{{ error checker
-
-class ErrorChecker(CSECachingMapperMixin, IdentityMapper):
- map_common_subexpression_uncached = \
- IdentityMapper.map_common_subexpression
-
- def __init__(self, mesh):
- self.mesh = mesh
-
- def map_operator_binding(self, expr):
- if isinstance(expr.op, op.DiffOperatorBase):
- if (self.mesh is not None
- and expr.op.xyz_axis >= self.mesh.ambient_dim):
- raise ValueError("optemplate tries to differentiate along a "
- "non-existent axis (e.g. Z in 2D)")
-
- # FIXME: Also check fluxes
- return IdentityMapper.map_operator_binding(self, expr)
-
- def map_normal(self, expr):
- if self.mesh is not None and expr.axis >= self.mesh.dimensions:
- raise ValueError("optemplate tries to differentiate along a "
- "non-existent axis (e.g. Z in 2D)")
-
- return expr
-
-# }}}
-
-
-# {{{ collectors for various symbolic operator components
-
-# To maintain deterministic output in code generation, these collectors return
-# OrderedSets. (As an example for why this is useful, the order of collected
-# values determines the names of intermediate variables. If variable names
-# aren't deterministic, loopy's kernel caching doesn't help us much across
-# runs.)
-
-class CollectorMixin(OperatorReducerMixin, LocalOpReducerMixin, FluxOpReducerMixin):
- def combine(self, values):
- from pytools import flatten
- return OrderedSet(flatten(values))
-
- def map_constant(self, expr, *args, **kwargs):
- return OrderedSet()
-
- map_variable = map_constant
- map_grudge_variable = map_constant
- map_function_symbol = map_constant
-
- map_ones = map_grudge_variable
- map_signed_face_ones = map_grudge_variable
- map_node_coordinate_component = map_grudge_variable
-
- map_operator = map_grudge_variable
-
-
-# I'm not sure this works still.
-#class GeometricFactorCollector(CollectorMixin, CombineMapper):
-# pass
-
-
-class BoundOperatorCollector(CSECachingMapperMixin, CollectorMixin, CombineMapper):
- def __init__(self, op_class):
- self.op_class = op_class
-
- map_common_subexpression_uncached = \
- CombineMapper.map_common_subexpression
-
- def map_operator_binding(self, expr):
- if isinstance(expr.op, self.op_class):
- result = OrderedSet([expr])
- else:
- result = OrderedSet()
-
- return result | CombineMapper.map_operator_binding(self, expr)
-
-
-class FluxExchangeCollector(CSECachingMapperMixin, CollectorMixin, CombineMapper):
- map_common_subexpression_uncached = \
- CombineMapper.map_common_subexpression
-
- def map_flux_exchange(self, expr):
- return OrderedSet([expr])
-
-# }}}
-
-
-# {{{ evaluation
-
-class Evaluator(pymbolic.mapper.evaluator.EvaluationMapper):
- pass
-
-
-class SymbolicEvaluator(pymbolic.mapper.evaluator.EvaluationMapper):
- def map_operator_binding(self, expr, *args, **kwargs):
- return expr.op(self.rec(expr.field, *args, **kwargs))
-
- def map_node_coordinate_component(self, expr, *args, **kwargs):
- return expr
-
- def map_call(self, expr, *args, **kwargs):
- return type(expr)(
- expr.function,
- tuple(self.rec(child, *args, **kwargs)
- for child in expr.parameters))
-
- def map_call_with_kwargs(self, expr, *args, **kwargs):
- return type(expr)(
- expr.function,
- tuple(self.rec(child, *args, **kwargs)
- for child in expr.parameters),
- {
- key: self.rec(val, *args, **kwargs)
- for key, val in expr.kw_parameters.items()}
- )
-
- def map_common_subexpression(self, expr):
- return type(expr)(self.rec(expr.child), expr.prefix, expr.scope)
-
-# }}}
-
-
-# vim: foldmethod=marker
diff --git a/grudge/symbolic/operators.py b/grudge/symbolic/operators.py
deleted file mode 100644
index 031ec60e..00000000
--- a/grudge/symbolic/operators.py
+++ /dev/null
@@ -1,862 +0,0 @@
-__copyright__ = "Copyright (C) 2008-2017 Andreas Kloeckner, Bogdan Enache"
-
-__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 sys import intern
-
-import numpy as np
-import pymbolic.primitives
-
-from typing import Tuple
-
-__doc__ = """
-
-Building blocks and mappers for operator expression trees.
-
-Basic Operators
-^^^^^^^^^^^^^^^
-
-.. autoclass:: Operator
-.. autoclass:: ElementwiseLinearOperator
-.. autoclass:: ProjectionOperator
-
-.. data:: project
-
-Reductions
-^^^^^^^^^^
-
-.. autoclass:: ElementwiseSumOperator
-.. autoclass:: ElementwiseMinOperator
-.. autoclass:: ElementwiseMaxOperator
-
-.. autoclass:: NodalReductionOperator
-.. autoclass:: NodalSum
-.. autoclass:: NodalMax
-.. autoclass:: NodalMin
-
-Differentiation
-^^^^^^^^^^^^^^^
-
-.. autoclass:: StrongFormDiffOperatorBase
-.. autoclass:: WeakFormDiffOperatorBase
-.. autoclass:: StiffnessOperator
-.. autoclass:: DiffOperator
-.. autoclass:: StiffnessTOperator
-.. autoclass:: MInvSTOperator
-
-.. autoclass:: RefDiffOperator
-.. autoclass:: RefStiffnessTOperator
-
-.. autofunction:: nabla
-.. autofunction:: minv_stiffness_t
-.. autofunction:: stiffness
-.. autofunction:: stiffness_t
-
-Mass Operators
-^^^^^^^^^^^^^^
-
-.. autoclass:: MassOperatorBase
-
-.. autoclass:: MassOperator
-.. autoclass:: InverseMassOperator
-.. autoclass:: FaceMassOperator
-
-.. autoclass:: RefMassOperator
-.. autoclass:: RefInverseMassOperator
-.. autoclass:: RefFaceMassOperator
-
-"""
-
-
-# {{{ base classes
-
-class Operator(pymbolic.primitives.Expression):
- """
- .. attribute:: dd_in
-
- an instance of :class:`~grudge.dof_desc.DOFDesc` describing the
- input discretization.
-
- .. attribute:: dd_out
-
- an instance of :class:`~grudge.dof_desc.DOFDesc` describing the
- output discretization.
- """
-
- def __init__(self, dd_in, dd_out):
- import grudge.dof_desc as dof_desc
- self.dd_in = dof_desc.as_dofdesc(dd_in)
- self.dd_out = dof_desc.as_dofdesc(dd_out)
-
- def stringifier(self):
- from grudge.symbolic.mappers import StringifyMapper
- return StringifyMapper
-
- def __call__(self, expr):
- from pytools.obj_array import obj_array_vectorize
- from grudge.tools import is_zero
-
- def bind_one(subexpr):
- if is_zero(subexpr):
- return subexpr
- else:
- from grudge.symbolic.primitives import OperatorBinding
- return OperatorBinding(self, subexpr)
-
- return obj_array_vectorize(bind_one, expr)
-
- def with_dd(self, dd_in=None, dd_out=None):
- """Return a copy of *self*, modified to the given DOF descriptors.
- """
- return type(self)(
- *self.__getinitargs__()[:-2],
- dd_in=dd_in or self.dd_in,
- dd_out=dd_out or self.dd_out)
-
- init_arg_names: Tuple[str, ...] = ("dd_in", "dd_out")
-
- def __getinitargs__(self):
- return (self.dd_in, self.dd_out,)
-
- def make_stringifier(self, originating_stringifier=None):
- from grudge.symbolic.mappers import StringifyMapper
- return StringifyMapper()
-# }}}
-
-
-class ElementwiseLinearOperator(Operator):
- def matrix(self, element_group):
- raise NotImplementedError
-
- mapper_method = intern("map_elementwise_linear")
-
-
-class ProjectionOperator(Operator):
- def __init__(self, dd_in, dd_out):
- super().__init__(dd_in, dd_out)
-
- def __call__(self, expr):
- from pytools.obj_array import obj_array_vectorize
-
- def project_one(subexpr):
- from pymbolic.primitives import is_constant
- if self.dd_in == self.dd_out:
- # no-op projection, go away
- return subexpr
- elif is_constant(subexpr):
- return subexpr
- else:
- from grudge.symbolic.primitives import OperatorBinding
- return OperatorBinding(self, subexpr)
-
- return obj_array_vectorize(project_one, expr)
-
- mapper_method = intern("map_projection")
-
-
-project = ProjectionOperator
-
-
-class InterpolationOperator(ProjectionOperator):
- def __init__(self, dd_in, dd_out):
- from warnings import warn
- warn("'InterpolationOperator' is deprecated, "
- "use 'ProjectionOperator' instead.",
- DeprecationWarning, stacklevel=2)
-
- super().__init__(dd_in, dd_out)
-
-
-def interp(dd_in, dd_out):
- from warnings import warn
- warn("using 'interp' is deprecated, use 'project' instead.",
- DeprecationWarning, stacklevel=2)
-
- return ProjectionOperator(dd_in, dd_out)
-
-
-# {{{ element reduction: sum, min, max
-
-class ElementwiseReductionOperator(Operator):
- def __init__(self, dd):
- super().__init__(dd_in=dd, dd_out=dd)
-
-
-class ElementwiseSumOperator(ElementwiseReductionOperator):
- """Returns a vector of DOFs with all entries on each element set
- to the sum of DOFs on that element.
- """
-
- mapper_method = intern("map_elementwise_sum")
-
-
-class ElementwiseMinOperator(ElementwiseReductionOperator):
- """Returns a vector of DOFs with all entries on each element set
- to the minimum of DOFs on that element.
- """
-
- mapper_method = intern("map_elementwise_min")
-
-
-class ElementwiseMaxOperator(ElementwiseReductionOperator):
- """Returns a vector of DOFs with all entries on each element set
- to the maximum of DOFs on that element.
- """
-
- mapper_method = intern("map_elementwise_max")
-
-# }}}
-
-
-# {{{ nodal reduction: sum, integral, max
-
-class NodalReductionOperator(Operator):
- def __init__(self, dd_in, dd_out=None):
- if dd_out is None:
- import grudge.dof_desc as dof_desc
- dd_out = dof_desc.DD_SCALAR
-
- assert dd_out.is_scalar()
-
- super().__init__(dd_out=dd_out, dd_in=dd_in)
-
-
-class NodalSum(NodalReductionOperator):
- mapper_method = intern("map_nodal_sum")
-
-
-class NodalMax(NodalReductionOperator):
- mapper_method = intern("map_nodal_max")
-
-
-class NodalMin(NodalReductionOperator):
- mapper_method = intern("map_nodal_min")
-
-# }}}
-
-
-# {{{ differentiation operators
-
-# {{{ global differentiation
-
-class DiffOperatorBase(Operator):
- def __init__(self, xyz_axis, dd_in=None, dd_out=None):
- import grudge.dof_desc as dof_desc
- if dd_in is None:
- dd_in = dof_desc.DD_VOLUME
-
- if dd_out is None:
- dd_out = dd_in.with_discr_tag(dof_desc.DISCR_TAG_BASE)
- else:
- dd_out = dof_desc.as_dofdesc(dd_out)
-
- if dd_out.uses_quadrature():
- raise ValueError("differentiation outputs are not on "
- "quadrature grids")
-
- super().__init__(dd_in, dd_out)
-
- self.xyz_axis = xyz_axis
-
- init_arg_names = ("xyz_axis", "dd_in", "dd_out")
-
- def __getinitargs__(self):
- return (self.xyz_axis, self.dd_in, self.dd_out)
-
- def equal_except_for_axis(self, other):
- return (type(self) == type(other)
- # first argument is always the axis
- and self.__getinitargs__()[1:] == other.__getinitargs__()[1:])
-
-
-class StrongFormDiffOperatorBase(DiffOperatorBase):
- pass
-
-
-class WeakFormDiffOperatorBase(DiffOperatorBase):
- pass
-
-
-class StiffnessOperator(StrongFormDiffOperatorBase):
- mapper_method = intern("map_stiffness")
-
-
-class DiffOperator(StrongFormDiffOperatorBase):
- mapper_method = intern("map_diff")
-
-
-class StiffnessTOperator(WeakFormDiffOperatorBase):
- mapper_method = intern("map_stiffness_t")
-
-
-class MInvSTOperator(WeakFormDiffOperatorBase):
- mapper_method = intern("map_minv_st")
-
-# }}}
-
-
-# {{{ reference-element differentiation
-
-class RefDiffOperatorBase(ElementwiseLinearOperator):
- def __init__(self, rst_axis, dd_in=None, dd_out=None):
- import grudge.dof_desc as dof_desc
- if dd_in is None:
- dd_in = dof_desc.DD_VOLUME
-
- if dd_out is None:
- dd_out = dd_in.with_discr_tag(dof_desc.DISCR_TAG_BASE)
-
- if dd_out.uses_quadrature():
- raise ValueError("differentiation outputs are not on "
- "quadrature grids")
-
- super().__init__(dd_in, dd_out)
-
- self.rst_axis = rst_axis
-
- init_arg_names = ("rst_axis", "dd_in", "dd_out")
-
- def __getinitargs__(self):
- return (self.rst_axis, self.dd_in, self.dd_out)
-
- def equal_except_for_axis(self, other):
- return (type(self) == type(other)
- # first argument is always the axis
- and self.__getinitargs__()[1:] == other.__getinitargs__()[1:])
-
-
-class RefDiffOperator(RefDiffOperatorBase):
- mapper_method = intern("map_ref_diff")
-
- @staticmethod
- def matrices(out_element_group, in_element_group):
- assert in_element_group == out_element_group
- from meshmode.discretization.poly_element import diff_matrices
- return diff_matrices(in_element_group)
-
-
-class RefStiffnessTOperator(RefDiffOperatorBase):
- mapper_method = intern("map_ref_stiffness_t")
-
- @staticmethod
- def matrices(out_elem_grp, in_elem_grp):
- if in_elem_grp == out_elem_grp:
- assert in_elem_grp.is_orthonormal_basis()
- from meshmode.discretization.poly_element import (mass_matrix,
- diff_matrices)
- mmat = mass_matrix(in_elem_grp)
- return [dmat.T.dot(mmat.T) for dmat in diff_matrices(in_elem_grp)]
-
- from modepy import vandermonde
- basis = out_elem_grp.basis_obj()
- vand = vandermonde(basis.functions, out_elem_grp.unit_nodes)
- grad_vand = vandermonde(basis.gradients, in_elem_grp.unit_nodes)
- vand_inv_t = np.linalg.inv(vand).T
-
- if not isinstance(grad_vand, tuple):
- # NOTE: special case for 1d
- grad_vand = (grad_vand,)
-
- weights = in_elem_grp.quadrature_rule().weights
- return np.einsum("c,bz,acz->abc", weights, vand_inv_t, grad_vand)
-
-
-# }}}
-
-# }}}
-
-
-# {{{ various elementwise linear operators
-
-class FilterOperator(ElementwiseLinearOperator):
- def __init__(self, mode_response_func, dd_in=None, dd_out=None):
- """
- :param mode_response_func: A function mapping
- ``(mode_tuple, local_discretization)`` to a float indicating the
- factor by which this mode is to be multiplied after filtering.
- (For example an instance of
- :class:`ExponentialFilterResponseFunction`.
- """
- import grudge.dof_desc as dof_desc
- if dd_in is None:
- dd_in = dof_desc.DD_VOLUME
-
- if dd_out is None:
- dd_out = dd_in
-
- if dd_in.uses_quadrature():
- raise ValueError("dd_in may not use quadrature")
- if dd_in != dd_out:
- raise ValueError("dd_in and dd_out must be identical")
-
- super().__init__(dd_in, dd_out)
-
- self.mode_response_func = mode_response_func
-
- def __getinitargs__(self):
- return (self.mode_response_func, self.dd_in, self.dd_out)
-
- def matrix(self, eg):
- # FIXME
- raise NotImplementedError()
-
- ldis = eg.local_discretization
-
- filter_coeffs = [self.mode_response_func(mid, ldis)
- for mid in ldis.generate_mode_identifiers()]
-
- # build filter matrix
- import numpy.linalg as la
- vdm = ldis.vandermonde()
- mat = np.asarray(
- la.solve(vdm.T, np.diag(filter_coeffs) @ vdm.T).T,
- order="C")
-
- return mat
-
-
-class AveragingOperator(ElementwiseLinearOperator):
- def __init__(self, dd_in=None, dd_out=None):
- import grudge.dof_desc as dof_desc
- if dd_in is None:
- dd_in = dof_desc.DD_VOLUME
-
- if dd_out is None:
- dd_out = dd_in
-
- if dd_in.uses_quadrature():
- raise ValueError("dd_in may not use quadrature")
- if dd_in != dd_out:
- raise ValueError("dd_in and dd_out must be identical")
-
- super().__init__(dd_in, dd_out)
-
- def matrix(self, eg):
- # average matrix, so that AVE*fields = cellaverage(fields)
- # see Hesthaven and Warburton page 227
-
- # FIXME
- raise NotImplementedError()
-
- mmat = eg.local_discretization.mass_matrix()
- standard_el_vol = np.sum(np.dot(mmat, np.ones(mmat.shape[0])))
- avg_mat_row = np.sum(mmat, 0)/standard_el_vol
-
- avg_mat = np.zeros((np.size(avg_mat_row), np.size(avg_mat_row)))
- avg_mat[:] = avg_mat_row
- return avg_mat
-
-
-class InverseVandermondeOperator(ElementwiseLinearOperator):
- def matrix(self, element_group):
- raise NotImplementedError() # FIXME
-
-
-class VandermondeOperator(ElementwiseLinearOperator):
- def matrix(self, element_group):
- raise NotImplementedError() # FIXME
-
-# }}}
-
-
-# {{{ mass operators
-
-class MassOperatorBase(Operator):
- """
- Inherits from :class:`Operator`.
-
- :attr:`~Operator.dd_in` and :attr:`~Operator.dd_out` may be surface or volume
- discretizations.
- """
-
- def __init__(self, dd_in=None, dd_out=None):
- import grudge.dof_desc as dof_desc
- if dd_in is None:
- dd_in = dof_desc.DD_VOLUME
- if dd_out is None:
- dd_out = dof_desc.DD_VOLUME
-
- super().__init__(dd_in, dd_out)
-
-
-class MassOperator(MassOperatorBase):
- mapper_method = intern("map_mass")
-
-
-class InverseMassOperator(MassOperatorBase):
- mapper_method = intern("map_inverse_mass")
-
-
-class RefMassOperatorBase(ElementwiseLinearOperator):
- pass
-
-
-class RefMassOperator(RefMassOperatorBase):
- @staticmethod
- def matrix(out_element_group, in_element_group):
- if out_element_group == in_element_group:
- from meshmode.discretization.poly_element import mass_matrix
- return mass_matrix(in_element_group)
-
- from modepy import vandermonde
- basis = out_element_group.basis_obj()
- vand = vandermonde(basis.functions, out_element_group.unit_nodes)
- o_vand = vandermonde(basis.functions, in_element_group.unit_nodes)
- vand_inv_t = np.linalg.inv(vand).T
-
- weights = in_element_group.quadrature_rule().weights
- return np.einsum("j,ik,jk->ij", weights, vand_inv_t, o_vand)
-
- mapper_method = intern("map_ref_mass")
-
-
-class RefInverseMassOperator(RefMassOperatorBase):
- @staticmethod
- def matrix(in_element_group, out_element_group):
- assert in_element_group == out_element_group
- import modepy as mp
- basis = in_element_group.basis_obj()
- return mp.inverse_mass_matrix(
- basis.functions,
- in_element_group.unit_nodes)
-
- mapper_method = intern("map_ref_inverse_mass")
-
-# }}}
-
-
-# {{{ boundary-related operators
-
-class OppositeInteriorFaceSwap(Operator):
- """
- .. attribute:: unique_id
-
- An integer identifying this specific instances of
- :class:`OppositePartitionFaceSwap` within an entire bound symbolic
- operator. Is assigned automatically by :func:`grudge.bind`
- if not already set by the user. This will become
- :class:`OppositePartitionFaceSwap.unique_id` in distributed
- runs.
- """
-
- def __init__(self, dd_in=None, dd_out=None, unique_id=None):
- from meshmode.discretization.connection import FACE_RESTR_INTERIOR
- import grudge.dof_desc as dof_desc
-
- if dd_in is None:
- dd_in = dof_desc.DOFDesc(FACE_RESTR_INTERIOR, None)
- if dd_out is None:
- dd_out = dd_in
-
- super().__init__(dd_in, dd_out)
- if self.dd_in.domain_tag is not FACE_RESTR_INTERIOR:
- raise ValueError("dd_in must be an interior faces domain")
- if self.dd_out != self.dd_in:
- raise ValueError("dd_out and dd_in must be identical")
-
- assert unique_id is None or isinstance(unique_id, int)
- self.unique_id = unique_id
-
- init_arg_names = ("dd_in", "dd_out", "unique_id")
-
- def __getinitargs__(self):
- return (self.dd_in, self.dd_out, self.unique_id)
-
- mapper_method = intern("map_opposite_interior_face_swap")
-
-
-class OppositePartitionFaceSwap(Operator):
- """
- .. attribute:: unique_id
-
- An integer corresponding to the :attr:`OppositeInteriorFaceSwap.unique_id`
- which led to the creation of this object. This integer is used as an
- MPI tag offset to keep different subexpressions apart in MPI traffic.
- """
- def __init__(self, dd_in=None, dd_out=None, unique_id=None):
- from meshmode.mesh import BTAG_PARTITION
- import grudge.dof_desc as dof_desc
-
- if dd_in is None and dd_out is None:
- raise ValueError("dd_in or dd_out must be specified")
- elif dd_in is None:
- dd_in = dd_out
- elif dd_out is None:
- dd_out = dd_in
-
- super().__init__(dd_in, dd_out)
- if not (isinstance(self.dd_in.domain_tag, dof_desc.DTAG_BOUNDARY)
- and isinstance(self.dd_in.domain_tag.tag, BTAG_PARTITION)):
- raise ValueError(
- "dd_in must be a partition boundary faces domain, not '%s'"
- % self.dd_in.domain_tag)
- if self.dd_out != self.dd_in:
- raise ValueError("dd_out and dd_in must be identical")
-
- self.i_remote_part = self.dd_in.domain_tag.tag.part_nr
-
- assert unique_id is None or isinstance(unique_id, int)
- self.unique_id = unique_id
-
- init_arg_names = ("dd_in", "dd_out", "unique_id")
-
- def __getinitargs__(self):
- return (self.dd_in, self.dd_out, self.unique_id)
-
- mapper_method = intern("map_opposite_partition_face_swap")
-
-
-class FaceMassOperatorBase(ElementwiseLinearOperator):
- def __init__(self, dd_in=None, dd_out=None):
- from meshmode.discretization.connection import FACE_RESTR_ALL
- import grudge.dof_desc as dof_desc
-
- if dd_in is None:
- dd_in = dof_desc.DOFDesc(FACE_RESTR_ALL, None)
-
- if dd_out is None or dd_out == "vol":
- dd_out = dof_desc.DOFDesc("vol", dof_desc.DISCR_TAG_BASE)
-
- if dd_out.uses_quadrature():
- raise ValueError("face mass operator outputs are not on "
- "quadrature grids")
-
- if not dd_out.is_volume():
- raise ValueError("dd_out must be a volume domain")
- if dd_in.domain_tag is not FACE_RESTR_ALL:
- raise ValueError("dd_in must be an interior faces domain")
-
- super().__init__(dd_in, dd_out)
-
-
-class FaceMassOperator(FaceMassOperatorBase):
- mapper_method = intern("map_face_mass_operator")
-
-
-class RefFaceMassOperator(ElementwiseLinearOperator):
- def matrix(self, afgrp, volgrp, dtype):
- nfaces = volgrp.mesh_el_group.nfaces
- assert afgrp.nelements == nfaces * volgrp.nelements
-
- matrix = np.empty(
- (volgrp.nunit_dofs,
- nfaces,
- afgrp.nunit_dofs),
- dtype=dtype)
-
- import modepy as mp
- from meshmode.discretization import ElementGroupWithBasis
- from meshmode.discretization.poly_element import \
- QuadratureSimplexElementGroup
-
- n = volgrp.order
- m = afgrp.order
- vol_basis = volgrp.basis_obj()
- faces = mp.faces_for_shape(volgrp.shape)
-
- for iface, face in enumerate(faces):
- # If the face group is defined on a higher-order
- # quadrature grid, use the underlying quadrature rule
- if isinstance(afgrp, QuadratureSimplexElementGroup):
- face_quadrature = afgrp.quadrature_rule()
- if face_quadrature.exact_to < m:
- raise ValueError(
- "The face quadrature rule is only exact for polynomials "
- f"of total degree {face_quadrature.exact_to}. Please "
- "ensure a quadrature rule is used that is at least "
- f"exact for degree {m}."
- )
- else:
- # NOTE: This handles the general case where
- # volume and surface quadrature rules may have different
- # integration orders
- face_quadrature = mp.quadrature_for_space(
- mp.space_for_shape(face, 2*max(n, m)),
- face
- )
-
- # If the group has a nodal basis and is unisolvent,
- # we use the basis on the face to compute the face mass matrix
- if (isinstance(afgrp, ElementGroupWithBasis)
- and afgrp.space.space_dim == afgrp.nunit_dofs):
-
- face_basis = afgrp.basis_obj()
-
- # Sanity check for face quadrature accuracy. Not integrating
- # degree N + M polynomials here is asking for a bad time.
- if face_quadrature.exact_to < m + n:
- raise ValueError(
- "The face quadrature rule is only exact for polynomials "
- f"of total degree {face_quadrature.exact_to}. Please "
- "ensure a quadrature rule is used that is at least "
- f"exact for degree {n+m}."
- )
-
- matrix[:, iface, :] = mp.nodal_mass_matrix_for_face(
- face, face_quadrature,
- face_basis.functions, vol_basis.functions,
- volgrp.unit_nodes, afgrp.unit_nodes,
- )
- else:
- # Otherwise, we use a routine that is purely quadrature-based
- # (no need for explicit face basis functions)
- matrix[:, iface, :] = mp.nodal_quad_mass_matrix_for_face(
- face,
- face_quadrature,
- vol_basis.functions,
- volgrp.unit_nodes,
- )
-
- return matrix
-
- mapper_method = intern("map_ref_face_mass_operator")
-
-# }}}
-
-
-# {{{ convenience functions for operator creation
-
-def nabla(dim):
- from pytools.obj_array import make_obj_array
- return make_obj_array(
- [DiffOperator(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, dd_in=None, dd_out=None):
- from pytools.obj_array import make_obj_array
- return make_obj_array(
- [StiffnessTOperator(i, dd_in, dd_out) for i in range(dim)])
-
-
-def integral(arg, dd=None):
- import grudge.symbolic.primitives as prim
- import grudge.dof_desc as dof_desc
-
- if dd is None:
- dd = dof_desc.DD_VOLUME
- dd = dof_desc.as_dofdesc(dd)
-
- return NodalSum(dd)(
- arg * prim.cse(
- MassOperator(dd_in=dd, dd_out=dd)(prim.Ones(dd)),
- "mass_quad_weights",
- prim.cse_scope.DISCRETIZATION))
-
-
-def norm(p, arg, dd=None):
- """
- :arg arg: is assumed to be a vector, i.e. have shape ``(n,)``.
- """
- import grudge.symbolic.primitives as prim
- import grudge.dof_desc as dof_desc
-
- if dd is None:
- dd = dof_desc.DD_VOLUME
- dd = dof_desc.as_dofdesc(dd)
-
- if p == 2:
- norm_squared = NodalSum(dd_in=dd)(
- arg * MassOperator()(arg))
-
- if isinstance(norm_squared, np.ndarray):
- if len(norm_squared.shape) != 1:
- raise NotImplementedError("can only take the norm of vectors")
-
- norm_squared = norm_squared.sum()
-
- return prim.sqrt(norm_squared)
-
- elif p == np.inf:
- result = NodalMax(dd_in=dd)(prim.fabs(arg))
-
- if isinstance(result, np.ndarray):
- if len(result.shape) != 1:
- raise NotImplementedError("can only take the norm of vectors")
-
- from pymbolic.primitives import Max
- result = Max(result)
-
- return result
-
- else:
- raise ValueError("unsupported value of p")
-
-
-def h_max_from_volume(ambient_dim, dim=None, dd=None):
- """Defines a characteristic length based on the volume of the elements.
- This length may not be representative if the elements have very high
- aspect ratios.
- """
-
- import grudge.symbolic.primitives as prim
- import grudge.dof_desc as dof_desc
-
- if dd is None:
- dd = dof_desc.DD_VOLUME
- dd = dof_desc.as_dofdesc(dd)
-
- if dim is None:
- dim = ambient_dim
-
- return NodalMax(dd_in=dd)(
- ElementwiseSumOperator(dd)(
- MassOperator(dd_in=dd, dd_out=dd)(prim.Ones(dd))
- )
- )**(1.0/dim)
-
-
-def h_min_from_volume(ambient_dim, dim=None, dd=None):
- """Defines a characteristic length based on the volume of the elements.
- This length may not be representative if the elements have very high
- aspect ratios.
- """
-
- import grudge.symbolic.primitives as prim
- import grudge.dof_desc as dof_desc
-
- if dd is None:
- dd = dof_desc.DD_VOLUME
- dd = dof_desc.as_dofdesc(dd)
-
- if dim is None:
- dim = ambient_dim
-
- return NodalMin(dd_in=dd)(
- ElementwiseSumOperator(dd)(
- MassOperator(dd_in=dd, dd_out=dd)(prim.Ones(dd))
- )
- )**(1.0/dim)
-
-# }}}
-
-# vim: foldmethod=marker
diff --git a/grudge/symbolic/primitives.py b/grudge/symbolic/primitives.py
deleted file mode 100644
index d95bc657..00000000
--- a/grudge/symbolic/primitives.py
+++ /dev/null
@@ -1,724 +0,0 @@
-"""Operator template language: primitives."""
-
-__copyright__ = "Copyright (C) 2008-2017 Andreas Kloeckner, Bogdan Enache"
-
-__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 sys import intern
-
-import numpy as np
-from pytools.obj_array import make_obj_array
-
-import pymbolic.primitives as prim
-from pymbolic.primitives import (
- Variable as VariableBase,
- CommonSubexpression as CommonSubexpressionBase,
- cse_scope as cse_scope_base,
- make_common_subexpression as cse)
-from pymbolic.geometric_algebra import MultiVector
-from grudge.trace_pair import TracePair
-
-
-class ExpressionBase(prim.Expression):
- def make_stringifier(self, originating_stringifier=None):
- from grudge.symbolic.mappers import StringifyMapper
- return StringifyMapper()
-
-
-__doc__ = """
-
-.. currentmodule:: grudge.sym
-
-Symbols
-^^^^^^^
-
-.. autoclass:: Variable
-.. autoclass:: ScalarVariable
-.. autoclass:: FunctionSymbol
-
-.. autofunction:: make_sym_array
-.. autofunction:: make_sym_mv
-
-.. function :: fabs(arg)
-.. function :: sqrt(arg)
-.. function :: exp(arg)
-.. function :: sin(arg)
-.. function :: cos(arg)
-.. function :: besesl_j(n, arg)
-.. function :: besesl_y(n, arg)
-
-Helpers
-^^^^^^^
-
-.. autoclass:: OperatorBinding
-.. autoclass:: PrioritizedSubexpression
-.. autoclass:: Ones
-
-Geometry data
-^^^^^^^^^^^^^
-
-.. autoclass:: NodeCoordinateComponent
-.. autofunction:: nodes
-.. autofunction:: mv_nodes
-.. autofunction:: forward_metric_derivative
-.. autofunction:: inverse_metric_derivative
-.. autofunction:: forward_metric_derivative_mat
-.. autofunction:: inverse_metric_derivative_mat
-.. autofunction:: first_fundamental_form
-.. autofunction:: inverse_first_fundamental_form
-.. autofunction:: inverse_surface_metric_derivative
-.. autofunction:: second_fundamental_form
-.. autofunction:: shape_operator
-.. autofunction:: pseudoscalar
-.. autofunction:: area_element
-.. autofunction:: mv_normal
-.. autofunction:: normal
-.. autofunction:: surface_normal
-.. autofunction:: summed_curvature
-.. autofunction:: mean_curvature
-
-Symbolic trace pair functions
-^^^^^^^^^^^^^^^^^^^^^^^^^^^^^
-
-.. autofunction:: int_tpair
-.. autofunction:: bv_tpair
-.. autofunction:: bdry_tpair
-"""
-
-
-# {{{ has-dof-desc mix-in
-
-class HasDOFDesc:
- """
- .. attribute:: dd
-
- an instance of :class:`grudge.dof_desc.DOFDesc` describing the
- discretization on which this property is given.
- """
-
- def __init__(self, *args, **kwargs):
- # The remaining arguments are passed to the chained superclass.
-
- if "dd" in kwargs:
- dd = kwargs.pop("dd")
- else:
- dd = args[-1]
- args = args[:-1]
-
- super().__init__(*args, **kwargs)
- self.dd = dd
-
- def __getinitargs__(self):
- return (self.dd,)
-
- def with_dd(self, dd):
- """Return a copy of *self*, modified to the given DOF descriptor.
- """
- return type(self)(*self.__getinitargs__())
-
-# }}}
-
-
-# {{{ variables
-
-class cse_scope(cse_scope_base): # noqa: N801
- DISCRETIZATION = "grudge_discretization"
-
-
-class Variable(HasDOFDesc, ExpressionBase, VariableBase):
- """A user-supplied input variable with a known
- :class:`grudge.dof_desc.DOFDesc`.
- """
- init_arg_names = ("name", "dd")
-
- def __init__(self, name, dd=None):
- import grudge.dof_desc as dof_desc
-
- if dd is None:
- dd = dof_desc.DD_VOLUME
-
- super().__init__(name, dd)
-
- def __getinitargs__(self):
- return (self.name, self.dd,)
-
- mapper_method = "map_grudge_variable"
-
-
-var = Variable
-
-
-class ScalarVariable(Variable):
- def __init__(self, name):
- import grudge.dof_desc as dof_desc
-
- super().__init__(name, dof_desc.DD_SCALAR)
-
-
-def make_sym_array(name, shape, dd=None):
- def var_factory(name):
- return Variable(name, dd)
-
- return prim.make_sym_array(name, shape, var_factory)
-
-
-def make_sym_mv(name, dim, var_factory=None):
- return MultiVector(
- make_sym_array(name, dim, var_factory))
-
-
-class FunctionSymbol(ExpressionBase, VariableBase):
- """A symbol to be used as the function argument of
- :class:`~pymbolic.primitives.Call`.
- """
-
- def __call__(self, *exprs):
- from pytools.obj_array import obj_array_vectorize_n_args
- return obj_array_vectorize_n_args(
- super().__call__, *exprs)
-
- mapper_method = "map_function_symbol"
-
-
-fabs = FunctionSymbol("fabs")
-sqrt = FunctionSymbol("sqrt")
-exp = FunctionSymbol("exp")
-sin = FunctionSymbol("sin")
-cos = FunctionSymbol("cos")
-atan2 = FunctionSymbol("atan2")
-bessel_j = FunctionSymbol("bessel_j")
-bessel_y = FunctionSymbol("bessel_y")
-
-# }}}
-
-
-# {{{ technical helpers
-
-class OperatorBinding(ExpressionBase):
- init_arg_names = ("op", "field")
-
- def __init__(self, op, field):
- self.op = op
- self.field = field
-
- mapper_method = intern("map_operator_binding")
-
- def __getinitargs__(self):
- return self.op, self.field
-
- def is_equal(self, other):
- return (other.__class__ == self.__class__
- and other.op == self.op
- and np.array_equal(other.field, self.field))
-
- def get_hash(self):
- from pytools.obj_array import obj_array_to_hashable
- return hash((self.__class__, self.op, obj_array_to_hashable(self.field)))
-
-
-class PrioritizedSubexpression(CommonSubexpressionBase):
- """When the optemplate-to-code transformation is performed,
- prioritized subexpressions work like common subexpression in
- that they are assigned their own separate identifier/register
- location. In addition to this behavior, prioritized subexpressions
- are evaluated with a settable priority, allowing the user to
- expedite or delay the evaluation of the subexpression.
- """
-
- def __init__(self, child, priority=0):
- super().__init__(child)
- self.priority = priority
-
- def __getinitargs__(self):
- return (self.child, self.priority)
-
- def get_extra_properties(self):
- return {"priority": self.priority}
-
-# }}}
-
-
-class Ones(HasDOFDesc, ExpressionBase):
- mapper_method = intern("map_ones")
-
-
-class _SignedFaceOnes(HasDOFDesc, ExpressionBase):
- """Produces DoFs on a face that are :math:`-1` if their corresponding
- face number is odd and :math:`+1` if it is even.
- *dd* must refer to a 0D (point-shaped) trace domain.
- This is based on the face order of
- :meth:`meshmode.mesh.MeshElementGroup.face_vertex_indices`.
-
- .. note::
-
- This is used as a hack to generate normals with the correct orientation
- in 1D problems, and so far only intended for this particular use cases.
- (If you can think of a better way, please speak up!)
- """
-
- def __init__(self, dd):
- import grudge.dof_desc as dof_desc
-
- dd = dof_desc.as_dofdesc(dd)
- assert dd.is_trace()
- super().__init__(dd)
-
- mapper_method = intern("map_signed_face_ones")
-
-
-# {{{ geometry data
-
-class DiscretizationProperty(HasDOFDesc, ExpressionBase):
- pass
-
-
-class NodeCoordinateComponent(DiscretizationProperty):
- def __init__(self, axis, dd=None):
- if not dd.is_discretized():
- raise ValueError("dd must be a discretization for "
- "NodeCoordinateComponent")
-
- super().__init__(dd)
- self.axis = axis
-
- assert dd.domain_tag is not None
-
- init_arg_names = ("axis", "dd")
-
- def __getinitargs__(self):
- return (self.axis, self.dd)
-
- mapper_method = intern("map_node_coordinate_component")
-
-
-def nodes(ambient_dim, dd=None):
- import grudge.dof_desc as dof_desc
-
- if dd is None:
- dd = dof_desc.DD_VOLUME
-
- dd = dof_desc.as_dofdesc(dd)
-
- return np.array([NodeCoordinateComponent(i, dd)
- for i in range(ambient_dim)], dtype=object)
-
-
-def mv_nodes(ambient_dim, dd=None):
- return MultiVector(nodes(ambient_dim, dd))
-
-
-def forward_metric_nth_derivative(xyz_axis, ref_axes, dd=None):
- r"""
- Pointwise metric derivatives representing repeated derivatives
-
- .. math::
-
- \frac{\partial^n x_{\mathrm{xyz\_axis}} }{\partial r_{\mathrm{ref\_axes}}}
-
- where *ref_axes* is a multi-index description.
-
- :arg ref_axes: a :class:`tuple` of tuples indicating indices of
- coordinate axes of the reference element to the number of derivatives
- which will be taken. For example, the value ``((0, 2), (1, 1))``
- indicates taking the second derivative with respect to the first
- axis and the first derivative with respect to the second
- axis. Each axis must occur only once and the tuple must be sorted
- by the axis index.
-
- May also be a singile integer *i*, which is viewed as equivalent
- to ``((i, 1),)``.
- """
- import grudge.dof_desc as dof_desc
-
- if isinstance(ref_axes, int):
- ref_axes = ((ref_axes, 1),)
-
- if not isinstance(ref_axes, tuple):
- raise ValueError("ref_axes must be a tuple")
-
- if tuple(sorted(ref_axes)) != ref_axes:
- raise ValueError("ref_axes must be sorted")
-
- if len(dict(ref_axes)) != len(ref_axes):
- raise ValueError("ref_axes must not contain an axis more than once")
-
- if dd is None:
- dd = dof_desc.DD_VOLUME
- inner_dd = dd.with_discr_tag(dof_desc.DISCR_TAG_BASE)
-
- from pytools import flatten
- flat_ref_axes = flatten([rst_axis] * n for rst_axis, n in ref_axes)
-
- from grudge.symbolic.operators import RefDiffOperator
- result = NodeCoordinateComponent(xyz_axis, inner_dd)
- for rst_axis in flat_ref_axes:
- result = RefDiffOperator(rst_axis, inner_dd)(result)
-
- if dd.uses_quadrature():
- from grudge.symbolic.operators import project
- result = project(inner_dd, dd)(result)
-
- return result
-
-
-def forward_metric_derivative(xyz_axis, rst_axis, dd=None):
- r"""
- Pointwise metric derivatives representing
-
- .. math::
-
- \frac{\partial x_{\mathrm{xyz\_axis}} }{\partial r_{\mathrm{rst\_axis}}}
- """
-
- return forward_metric_nth_derivative(xyz_axis, rst_axis, dd=dd)
-
-
-def forward_metric_derivative_vector(ambient_dim, rst_axis, dd=None):
- return make_obj_array([
- forward_metric_derivative(i, rst_axis, dd=dd)
- for i in range(ambient_dim)])
-
-
-def forward_metric_derivative_mv(ambient_dim, rst_axis, dd=None):
- return MultiVector(
- forward_metric_derivative_vector(ambient_dim, rst_axis, dd=dd))
-
-
-def parametrization_derivative(ambient_dim, dim=None, dd=None):
- if dim is None:
- dim = ambient_dim
-
- if dim == 0:
- from pymbolic.geometric_algebra import get_euclidean_space
- return MultiVector(_SignedFaceOnes(dd),
- space=get_euclidean_space(ambient_dim))
-
- from pytools import product
- return product(
- forward_metric_derivative_mv(ambient_dim, rst_axis, dd)
- for rst_axis in range(dim))
-
-
-def inverse_metric_derivative(rst_axis, xyz_axis, ambient_dim, dim=None,
- dd=None):
- if dim is None:
- dim = ambient_dim
-
- if dim != ambient_dim:
- raise ValueError("not clear what inverse_metric_derivative means if "
- "the derivative matrix is not square")
-
- par_vecs = [
- forward_metric_derivative_mv(ambient_dim, rst, dd)
- for rst in range(dim)]
-
- # Yay Cramer's rule! (o rly?)
- from functools import reduce, partial
- from operator import xor as outerprod_op
- outerprod = partial(reduce, outerprod_op)
-
- def outprod_with_unit(i, at):
- unit_vec = np.zeros(dim)
- unit_vec[i] = 1
-
- vecs = par_vecs[:]
- vecs[at] = MultiVector(unit_vec)
-
- return outerprod(vecs)
-
- volume_pseudoscalar_inv = cse(outerprod(
- forward_metric_derivative_mv(
- ambient_dim, rst_axis, dd=dd)
- for rst_axis in range(dim)).inv())
-
- return cse(
- (outprod_with_unit(xyz_axis, rst_axis)
- * volume_pseudoscalar_inv).as_scalar(),
- prefix="dr%d_dx%d" % (rst_axis, xyz_axis),
- scope=cse_scope.DISCRETIZATION)
-
-
-def forward_metric_derivative_mat(ambient_dim, dim=None, dd=None):
- if dim is None:
- dim = ambient_dim
-
- result = np.zeros((ambient_dim, dim), dtype=object)
- for j in range(dim):
- result[:, j] = forward_metric_derivative_vector(ambient_dim, j, dd=dd)
-
- return result
-
-
-def inverse_metric_derivative_mat(ambient_dim, dim=None, dd=None):
- if dim is None:
- dim = ambient_dim
-
- result = np.zeros((ambient_dim, dim), dtype=object)
- for i in range(dim):
- for j in range(ambient_dim):
- result[i, j] = inverse_metric_derivative(
- i, j, ambient_dim, dim, dd=dd)
-
- return result
-
-
-def first_fundamental_form(ambient_dim, dim=None, dd=None):
- if dim is None:
- dim = ambient_dim - 1
-
- mder = forward_metric_derivative_mat(ambient_dim, dim=dim, dd=dd)
- return cse(mder.T.dot(mder), "form1_mat", cse_scope.DISCRETIZATION)
-
-
-def inverse_first_fundamental_form(ambient_dim, dim=None, dd=None):
- if dim is None:
- dim = ambient_dim - 1
-
- if ambient_dim == dim:
- inv_mder = inverse_metric_derivative_mat(ambient_dim, dim=dim, dd=dd)
- inv_form1 = inv_mder.dot(inv_mder.T)
- else:
- form1 = first_fundamental_form(ambient_dim, dim, dd)
-
- if dim == 1:
- inv_form1 = np.array([[1.0/form1[0, 0]]], dtype=object)
- elif dim == 2:
- (E, F), (_, G) = form1 # noqa: N806
- inv_form1 = 1.0 / (E * G - F * F) * np.array([
- [G, -F], [-F, E]
- ], dtype=object)
- else:
- raise ValueError("%dD surfaces not supported" % dim)
-
- return cse(inv_form1, "inv_form1_mat", cse_scope.DISCRETIZATION)
-
-
-def inverse_surface_metric_derivative(rst_axis, xyz_axis,
- ambient_dim, dim=None, dd=None):
- if dim is None:
- dim = ambient_dim - 1
-
- if ambient_dim == dim:
- return inverse_metric_derivative(rst_axis, xyz_axis,
- ambient_dim, dim=dim, dd=dd)
- else:
- inv_form1 = inverse_first_fundamental_form(ambient_dim, dim=dim, dd=dd)
- imd = sum(
- inv_form1[rst_axis, d]*forward_metric_derivative(xyz_axis, d, dd=dd)
- for d in range(dim))
-
- return cse(imd,
- prefix="ds%d_dx%d" % (rst_axis, xyz_axis),
- scope=cse_scope.DISCRETIZATION)
-
-
-def second_fundamental_form(ambient_dim, dim=None, dd=None):
- if dim is None:
- dim = ambient_dim - 1
-
- normal = surface_normal(ambient_dim, dim=dim, dd=dd).as_vector()
- if dim == 1:
- second_ref_axes = [((0, 2),)]
- elif dim == 2:
- second_ref_axes = [((0, 2),), ((0, 1), (1, 1)), ((1, 2),)]
- else:
- raise ValueError("%dD surfaces not supported" % dim)
-
- from pytools import flatten
- form2 = np.empty((dim, dim), dtype=object)
- for ref_axes in second_ref_axes:
- i, j = flatten([rst_axis] * n for rst_axis, n in ref_axes)
-
- ruv = np.array([
- forward_metric_nth_derivative(xyz_axis, ref_axes, dd=dd)
- for xyz_axis in range(ambient_dim)])
- form2[i, j] = form2[j, i] = normal.dot(ruv)
-
- return cse(form2, "form2_mat", cse_scope.DISCRETIZATION)
-
-
-def shape_operator(ambient_dim, dim=None, dd=None):
- if dim is None:
- dim = ambient_dim - 1
-
- inv_form1 = inverse_first_fundamental_form(ambient_dim, dim=dim, dd=dd)
- form2 = second_fundamental_form(ambient_dim, dim=dim, dd=dd)
-
- return cse(-form2.dot(inv_form1), "shape_operator", cse_scope.DISCRETIZATION)
-
-
-def pseudoscalar(ambient_dim, dim=None, dd=None):
- if dim is None:
- dim = ambient_dim
-
- return parametrization_derivative(ambient_dim, dim, dd=dd).project_max_grade()
-
-
-def area_element(ambient_dim, dim=None, dd=None):
- return cse(
- sqrt(
- pseudoscalar(ambient_dim, dim, dd=dd)
- .norm_squared()),
- "area_element", cse_scope.DISCRETIZATION)
-
-
-def surface_normal(ambient_dim, dim=None, dd=None):
- import grudge.dof_desc as dof_desc
-
- dd = dof_desc.as_dofdesc(dd)
- if dim is None:
- dim = ambient_dim - 1
-
- # NOTE: Don't be tempted to add a sign here. As it is, it produces
- # exterior normals for positively oriented curves.
-
- pder = pseudoscalar(ambient_dim, dim, dd=dd) \
- / area_element(ambient_dim, dim, dd=dd)
-
- # Dorst Section 3.7.2
- return cse(pder << pder.I.inv(),
- "surface_normal",
- cse_scope.DISCRETIZATION)
-
-
-def mv_normal(dd, ambient_dim, dim=None):
- """Exterior unit normal as a :class:`~pymbolic.geometric_algebra.MultiVector`."""
- import grudge.dof_desc as dof_desc
-
- dd = dof_desc.as_dofdesc(dd)
- if not dd.is_trace():
- raise ValueError("may only request normals on boundaries")
-
- if dim is None:
- dim = ambient_dim - 1
-
- if dim == ambient_dim - 1:
- return surface_normal(ambient_dim, dim=dim, dd=dd)
-
- # NOTE: In the case of (d - 2)-dimensional curves, we don't really have
- # enough information on the face to decide what an "exterior face normal"
- # is (e.g the "normal" to a 1D curve in 3D space is actually a
- # "normal plane")
- #
- # The trick done here is that we take the surface normal, move it to the
- # face and then take a cross product with the face tangent to get the
- # correct exterior face normal vector.
- assert dim == ambient_dim - 2
-
- from grudge.symbolic.operators import project
- import grudge.dof_desc as dof_desc
-
- volm_normal = (
- surface_normal(ambient_dim, dim=dim + 1,
- dd=dof_desc.DD_VOLUME)
- .map(project(dof_desc.DD_VOLUME, dd)))
- pder = pseudoscalar(ambient_dim, dim, dd=dd)
-
- mv = cse(-(volm_normal ^ pder) << volm_normal.I.inv(),
- "face_normal",
- cse_scope.DISCRETIZATION)
-
- return cse(mv / sqrt(mv.norm_squared()),
- "unit_face_normal",
- cse_scope.DISCRETIZATION)
-
-
-def normal(dd, ambient_dim, dim=None):
- return mv_normal(dd, ambient_dim, dim).as_vector()
-
-
-def summed_curvature(ambient_dim, dim=None, dd=None):
- """Sum of the principal curvatures"""
-
- if dim is None:
- dim = ambient_dim - 1
-
- if ambient_dim == 1:
- return 0.0
-
- if ambient_dim == dim:
- return 0.0
-
- op = shape_operator(ambient_dim, dim=dim, dd=dd)
- return cse(np.trace(op), "summed_curvature", cse_scope.DISCRETIZATION)
-
-
-def mean_curvature(ambient_dim, dim=None, dd=None):
- """Averaged (by dimension) sum of the principal curvatures."""
- return 1.0 / (ambient_dim-1.0) * summed_curvature(ambient_dim, dim=dim, dd=dd)
-
-# }}}
-
-
-# {{{ Symbolic trace pair functions
-
-def int_tpair(expression, qtag=None, from_dd=None):
- from meshmode.discretization.connection import FACE_RESTR_INTERIOR
- from grudge.symbolic.operators import project, OppositeInteriorFaceSwap
- import grudge.dof_desc as dof_desc
-
- if from_dd is None:
- from_dd = dof_desc.DD_VOLUME
- assert not from_dd.uses_quadrature()
-
- trace_dd = dof_desc.DOFDesc(FACE_RESTR_INTERIOR, qtag)
- if from_dd.domain_tag == trace_dd.domain_tag:
- i = expression
- else:
- i = project(from_dd, trace_dd.with_qtag(None))(expression)
- e = cse(OppositeInteriorFaceSwap()(i))
-
- if trace_dd.uses_quadrature():
- i = cse(project(trace_dd.with_qtag(None), trace_dd)(i))
- e = cse(project(trace_dd.with_qtag(None), trace_dd)(e))
-
- return TracePair(trace_dd, interior=i, exterior=e)
-
-
-def bdry_tpair(dd, interior, exterior):
- """
- :arg interior: an expression that already lives on the boundary
- representing the interior value to be used
- for the flux.
- :arg exterior: an expression that already lives on the boundary
- representing the exterior value to be used
- for the flux.
- """
- return TracePair(dd, interior=interior, exterior=exterior)
-
-
-def bv_tpair(dd, interior, exterior):
- """
- :arg interior: an expression that lives in the volume
- and will be restricted to the boundary identified by
- *tag* before use.
- :arg exterior: an expression that already lives on the boundary
- representing the exterior value to be used
- for the flux.
- """
- from grudge.symbolic.operators import project
- interior = cse(project("vol", dd)(interior))
- return TracePair(dd, interior=interior, exterior=exterior)
-
-# }}}
-
-
-# vim: foldmethod=marker
diff --git a/grudge/symbolic/tools.py b/grudge/symbolic/tools.py
deleted file mode 100644
index 61776f1b..00000000
--- a/grudge/symbolic/tools.py
+++ /dev/null
@@ -1,141 +0,0 @@
-"""Operator templates: extra bits of functionality."""
-
-__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 numpy as np
-
-
-# {{{ symbolic operator tools
-
-def is_scalar(expr):
- from grudge import sym
- from grudge.dof_desc import DTAG_SCALAR
-
- if isinstance(expr, sym.Variable) and expr.dd.domain_tag is DTAG_SCALAR:
- return True
-
- return isinstance(expr, (int, float, complex))
-
-
-def split_sym_operator_for_multirate(state_vector, sym_operator,
- index_groups):
- class IndexGroupKillerSubstMap:
- def __init__(self, kill_set):
- self.kill_set = kill_set
-
- def __call__(self, expr):
- if expr in kill_set:
- return 0
- else:
- return None
-
- # make IndexGroupKillerSubstMap that kill everything
- # *except* what's in that index group
- killers = []
- for i in range(len(index_groups)):
- kill_set = set()
- for j in range(len(index_groups)):
- if i != j:
- kill_set |= set(index_groups[j])
-
- killers.append(IndexGroupKillerSubstMap(kill_set))
-
- from grudge.symbolic.mappers import \
- SubstitutionMapper, \
- CommutativeConstantFoldingMapper
-
- return [
- CommutativeConstantFoldingMapper()(
- SubstitutionMapper(killer)(
- sym_operator[ig]))
- for ig in index_groups
- for killer in killers]
-
-
-def ptwise_mul(a, b):
- from pytools.obj_array import log_shape
- a_log_shape = log_shape(a)
- b_log_shape = log_shape(b)
-
- from pytools import indices_in_shape
-
- if a_log_shape == ():
- result = np.empty(b_log_shape, dtype=object)
- for i in indices_in_shape(b_log_shape):
- result[i] = a*b[i]
- elif b_log_shape == ():
- result = np.empty(a_log_shape, dtype=object)
- for i in indices_in_shape(a_log_shape):
- result[i] = a[i]*b
- else:
- raise ValueError("ptwise_mul can't handle two non-scalars")
-
- return result
-
-
-def ptwise_dot(logdims1, logdims2, a1, a2):
- a1_log_shape = a1.shape[:logdims1]
- a2_log_shape = a1.shape[:logdims2]
-
- assert a1_log_shape[-1] == a2_log_shape[0]
- len_k = a2_log_shape[0]
-
- result = np.empty(a1_log_shape[:-1]+a2_log_shape[1:], dtype=object)
-
- from pytools import indices_in_shape
- for a1_i in indices_in_shape(a1_log_shape[:-1]):
- for a2_i in indices_in_shape(a2_log_shape[1:]):
- result[a1_i+a2_i] = sum(
- a1[a1_i+(k,)] * a2[(k,)+a2_i]
- for k in range(len_k)
- )
-
- if result.shape == ():
- return result[()]
- else:
- return result
-
-# }}}
-
-
-# {{{ pretty printing
-
-def pretty(sym_operator):
- from grudge.symbolic.mappers import PrettyStringifyMapper
-
- stringify_mapper = PrettyStringifyMapper()
- from pymbolic.mapper.stringifier import PREC_NONE
- result = stringify_mapper(sym_operator, PREC_NONE)
-
- splitter = "="*75 + "\n"
-
- cse_strs = stringify_mapper.get_cse_strings()
- if cse_strs:
- result = "\n".join(cse_strs)+"\n"+splitter+result
-
- return result
-
-# }}}
-
-
-# vim: foldmethod=marker
diff --git a/grudge/trace_pair.py b/grudge/trace_pair.py
index f84dfdc5..1fd5404b 100644
--- a/grudge/trace_pair.py
+++ b/grudge/trace_pair.py
@@ -105,11 +105,6 @@ class TracePair:
.. automethod:: __getattr__
.. automethod:: __getitem__
.. automethod:: __len__
-
- .. note::
-
- :class:`TracePair` is currently used both by the symbolic (deprecated)
- and the current interfaces, with symbolic information or concrete data.
"""
dd: dof_desc.DOFDesc
diff --git a/test/test_grudge.py b/test/test_grudge.py
index 7e469d46..6a42167a 100644
--- a/test/test_grudge.py
+++ b/test/test_grudge.py
@@ -946,6 +946,7 @@ def test_improvement_quadrature(actx_factory, order):
# {{{ bessel
+@pytest.mark.xfail
def test_bessel(actx_factory):
actx = actx_factory()
@@ -964,7 +965,7 @@ def test_bessel(actx_factory):
# FIXME: Bessel functions need to brought out of the symbolic
# layer. Related issue: https://github.com/inducer/grudge/issues/93
def bessel_j(actx, n, r):
- from grudge import sym, bind
+ from grudge import sym, bind # pylint: disable=no-name-in-module
return bind(dcoll, sym.bessel_j(n, sym.var("r")))(actx, r=r)
# https://dlmf.nist.gov/10.6.1
diff --git a/test/test_grudge_sym_old.py b/test/test_grudge_sym_old.py
deleted file mode 100644
index f211a316..00000000
--- a/test/test_grudge_sym_old.py
+++ /dev/null
@@ -1,943 +0,0 @@
-__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
-
-from meshmode.dof_array import flatten, thaw
-import meshmode.mesh.generation as mgen
-
-from pytools.obj_array import flat_obj_array, make_obj_array
-
-from grudge import sym, bind, DiscretizationCollection
-
-import grudge.dof_desc as dof_desc
-
-import pytest
-
-from grudge.array_context import PytestPyOpenCLArrayContextFactoryWithHostScalars
-from arraycontext import pytest_generate_tests_for_array_contexts
-pytest_generate_tests = pytest_generate_tests_for_array_contexts(
- [PytestPyOpenCLArrayContextFactoryWithHostScalars])
-
-import logging
-
-logger = logging.getLogger(__name__)
-
-
-# {{{ inverse metric
-
-@pytest.mark.parametrize("dim", [2, 3])
-def test_inverse_metric(actx_factory, dim):
- actx = actx_factory()
-
- mesh = mgen.generate_regular_rect_mesh(a=(-0.5,)*dim, b=(0.5,)*dim,
- nelements_per_axis=(6,)*dim, order=4)
-
- def m(x):
- result = np.empty_like(x)
- result[0] = (
- 1.5*x[0] + np.cos(x[0])
- + 0.1*np.sin(10*x[1]))
- result[1] = (
- 0.05*np.cos(10*x[0])
- + 1.3*x[1] + np.sin(x[1]))
- if len(x) == 3:
- result[2] = x[2]
- return result
-
- from meshmode.mesh.processing import map_mesh
- mesh = map_mesh(mesh, m)
-
- discr = DiscretizationCollection(actx, mesh, order=4)
-
- sym_op = (
- sym.forward_metric_derivative_mat(mesh.dim)
- .dot(
- sym.inverse_metric_derivative_mat(mesh.dim)
- )
- .reshape(-1))
-
- op = bind(discr, sym_op)
- mat = op(actx).reshape(mesh.dim, mesh.dim)
-
- for i in range(mesh.dim):
- for j in range(mesh.dim):
- tgt = 1 if i == j else 0
-
- err = actx.np.linalg.norm(mat[i, j] - tgt, ord=np.inf)
- logger.info("error[%d, %d]: %.5e", i, j, err)
- assert err < 1.0e-12, (i, j, err)
-
-# }}}
-
-
-# {{{ mass operator trig integration
-
-@pytest.mark.parametrize("ambient_dim", [1, 2, 3])
-@pytest.mark.parametrize("discr_tag", [dof_desc.DISCR_TAG_BASE,
- dof_desc.DISCR_TAG_QUAD])
-def test_mass_mat_trig(actx_factory, ambient_dim, discr_tag):
- """Check the integral of some trig functions on an interval using the mass
- matrix.
- """
- actx = actx_factory()
-
- nel_1d = 16
- order = 4
-
- a = -4.0 * np.pi
- b = +9.0 * np.pi
- true_integral = 13*np.pi/2 * (b - a)**(ambient_dim - 1)
-
- from meshmode.discretization.poly_element import QuadratureSimplexGroupFactory
- dd_quad = dof_desc.DOFDesc(dof_desc.DTAG_VOLUME_ALL, discr_tag)
- if discr_tag is dof_desc.DISCR_TAG_BASE:
- discr_tag_to_group_factory = {}
- else:
- discr_tag_to_group_factory = {
- discr_tag: QuadratureSimplexGroupFactory(order=2*order)
- }
-
- mesh = mgen.generate_regular_rect_mesh(
- a=(a,)*ambient_dim, b=(b,)*ambient_dim,
- nelements_per_axis=(nel_1d,)*ambient_dim, order=1)
- discr = DiscretizationCollection(
- actx, mesh, order=order,
- discr_tag_to_group_factory=discr_tag_to_group_factory
- )
-
- def _get_variables_on(dd):
- sym_f = sym.var("f", dd=dd)
- sym_x = sym.nodes(ambient_dim, dd=dd)
- sym_ones = sym.Ones(dd)
-
- return sym_f, sym_x, sym_ones
-
- sym_f, sym_x, sym_ones = _get_variables_on(dof_desc.DD_VOLUME)
- f_volm = actx.to_numpy(flatten(bind(discr, sym.cos(sym_x[0])**2)(actx)))
- ones_volm = actx.to_numpy(flatten(bind(discr, sym_ones)(actx)))
-
- sym_f, sym_x, sym_ones = _get_variables_on(dd_quad)
- f_quad = bind(discr, sym.cos(sym_x[0])**2)(actx)
- ones_quad = bind(discr, sym_ones)(actx)
-
- mass_op = bind(discr, sym.MassOperator(dd_quad, dof_desc.DD_VOLUME)(sym_f))
-
- num_integral_1 = np.dot(ones_volm, actx.to_numpy(flatten(mass_op(f=f_quad))))
- err_1 = abs(num_integral_1 - true_integral)
- assert err_1 < 9e-9, err_1
-
- num_integral_2 = np.dot(f_volm, actx.to_numpy(flatten(mass_op(f=ones_quad))))
- err_2 = abs(num_integral_2 - true_integral)
- assert err_2 < 5e-9, err_2
-
- if discr_tag is dof_desc.DISCR_TAG_BASE:
- # NOTE: `integral` always makes a square mass matrix and
- # `QuadratureSimplexGroupFactory` does not have a `mass_matrix` method.
- num_integral_3 = bind(discr,
- sym.integral(sym_f, dd=dd_quad))(f=f_quad)
- err_3 = abs(num_integral_3 - true_integral)
- assert err_3 < 5.0e-10, err_3
-
-# }}}
-
-
-# {{{ mass operator surface area
-
-def _ellipse_surface_area(radius, aspect_ratio):
- # https://docs.scipy.org/doc/scipy/reference/generated/scipy.special.ellipe.html
- eccentricity = 1.0 - (1/aspect_ratio)**2
-
- if abs(aspect_ratio - 2.0) < 1.0e-14:
- # NOTE: hardcoded value so we don't need scipy for the test
- ellip_e = 1.2110560275684594
- else:
- from scipy.special import ellipe # pylint: disable=no-name-in-module
- ellip_e = ellipe(eccentricity)
-
- return 4.0 * radius * ellip_e
-
-
-def _spheroid_surface_area(radius, aspect_ratio):
- # https://en.wikipedia.org/wiki/Ellipsoid#Surface_area
- a = 1.0
- c = aspect_ratio
-
- if a < c:
- e = np.sqrt(1.0 - (a/c)**2)
- return 2.0 * np.pi * radius**2 * (1.0 + (c/a) / e * np.arcsin(e))
- else:
- e = np.sqrt(1.0 - (c/a)**2)
- return 2.0 * np.pi * radius**2 * (1 + (c/a)**2 / e * np.arctanh(e))
-
-
-@pytest.mark.parametrize("name", [
- "2-1-ellipse", "spheroid", "box2d", "box3d"
- ])
-def test_mass_surface_area(actx_factory, name):
- actx = actx_factory()
-
- # {{{ cases
-
- if name == "2-1-ellipse":
- from mesh_data import EllipseMeshBuilder
- builder = EllipseMeshBuilder(radius=3.1, aspect_ratio=2.0)
- surface_area = _ellipse_surface_area(builder.radius, builder.aspect_ratio)
- elif name == "spheroid":
- from mesh_data import SpheroidMeshBuilder
- builder = SpheroidMeshBuilder()
- surface_area = _spheroid_surface_area(builder.radius, builder.aspect_ratio)
- elif name == "box2d":
- from mesh_data import BoxMeshBuilder
- builder = BoxMeshBuilder(ambient_dim=2)
- surface_area = 1.0
- elif name == "box3d":
- from mesh_data import BoxMeshBuilder
- builder = BoxMeshBuilder(ambient_dim=3)
- surface_area = 1.0
- else:
- raise ValueError("unknown geometry name: %s" % name)
-
- # }}}
-
- # {{{ convergence
-
- from pytools.convergence import EOCRecorder
- eoc = EOCRecorder()
-
- for resolution in builder.resolutions:
- mesh = builder.get_mesh(resolution, builder.mesh_order)
- discr = DiscretizationCollection(actx, mesh, order=builder.order)
- volume_discr = discr.discr_from_dd(dof_desc.DD_VOLUME)
-
- logger.info("ndofs: %d", volume_discr.ndofs)
- logger.info("nelements: %d", volume_discr.mesh.nelements)
-
- # {{{ compute surface area
-
- dd = dof_desc.DD_VOLUME
- sym_op = sym.NodalSum(dd)(sym.MassOperator(dd, dd)(sym.Ones(dd)))
- approx_surface_area = bind(discr, sym_op)(actx)
-
- logger.info("surface: got {:.5e} / expected {:.5e}".format(
- approx_surface_area, surface_area))
- area_error = abs(approx_surface_area - surface_area) / abs(surface_area)
-
- # }}}
-
- h_max = bind(discr, sym.h_max_from_volume(
- discr.ambient_dim, dim=discr.dim, dd=dd))(actx)
- eoc.add_data_point(h_max, area_error + 1.0e-16)
-
- # }}}
-
- logger.info("surface area error\n%s", str(eoc))
-
- assert eoc.max_error() < 1.0e-14 \
- or eoc.order_estimate() > builder.order
-
-# }}}
-
-
-# {{{ surface mass inverse
-
-@pytest.mark.parametrize("name", ["2-1-ellipse", "spheroid"])
-def test_surface_mass_operator_inverse(actx_factory, name):
- actx = actx_factory()
-
- # {{{ cases
-
- if name == "2-1-ellipse":
- from mesh_data import EllipseMeshBuilder
- builder = EllipseMeshBuilder(radius=3.1, aspect_ratio=2.0)
- elif name == "spheroid":
- from mesh_data import SpheroidMeshBuilder
- builder = SpheroidMeshBuilder()
- else:
- raise ValueError("unknown geometry name: %s" % name)
-
- # }}}
-
- # {{{ convergence
-
- from pytools.convergence import EOCRecorder
- eoc = EOCRecorder()
-
- for resolution in builder.resolutions:
- mesh = builder.get_mesh(resolution, builder.mesh_order)
- discr = DiscretizationCollection(actx, mesh, order=builder.order)
- volume_discr = discr.discr_from_dd(dof_desc.DD_VOLUME)
-
- logger.info("ndofs: %d", volume_discr.ndofs)
- logger.info("nelements: %d", volume_discr.mesh.nelements)
-
- # {{{ compute inverse mass
-
- dd = dof_desc.DD_VOLUME
- sym_f = sym.cos(4.0 * sym.nodes(mesh.ambient_dim, dd)[0])
- sym_op = sym.InverseMassOperator(dd, dd)(
- sym.MassOperator(dd, dd)(sym.var("f")))
-
- f = bind(discr, sym_f)(actx)
- f_inv = bind(discr, sym_op)(actx, f=f)
-
- inv_error = bind(discr,
- sym.norm(2, sym.var("x") - sym.var("y"))
- / sym.norm(2, sym.var("y")))(actx, x=f_inv, y=f)
-
- # }}}
-
- h_max = bind(discr, sym.h_max_from_volume(
- discr.ambient_dim, dim=discr.dim, dd=dd))(actx)
- eoc.add_data_point(h_max, inv_error)
-
- # }}}
-
- logger.info("inverse mass error\n%s", str(eoc))
-
- # NOTE: both cases give 1.0e-16-ish at the moment, but just to be on the
- # safe side, choose a slightly larger tolerance
- assert eoc.max_error() < 1.0e-14
-
-# }}}
-
-
-# {{{ surface face normal orthogonality
-
-@pytest.mark.parametrize("mesh_name", ["2-1-ellipse", "spheroid"])
-def test_face_normal_surface(actx_factory, mesh_name):
- """Check that face normals are orthogonal to the surface normal"""
- actx = actx_factory()
-
- # {{{ geometry
-
- if mesh_name == "2-1-ellipse":
- from mesh_data import EllipseMeshBuilder
- builder = EllipseMeshBuilder(radius=3.1, aspect_ratio=2.0)
- elif mesh_name == "spheroid":
- from mesh_data import SpheroidMeshBuilder
- builder = SpheroidMeshBuilder()
- else:
- raise ValueError("unknown mesh name: %s" % mesh_name)
-
- mesh = builder.get_mesh(builder.resolutions[0], builder.mesh_order)
- discr = DiscretizationCollection(actx, mesh, order=builder.order)
-
- volume_discr = discr.discr_from_dd(dof_desc.DD_VOLUME)
- logger.info("ndofs: %d", volume_discr.ndofs)
- logger.info("nelements: %d", volume_discr.mesh.nelements)
-
- # }}}
-
- # {{{ symbolic
- from meshmode.discretization.connection import FACE_RESTR_INTERIOR
-
- dv = dof_desc.DD_VOLUME
- df = dof_desc.as_dofdesc(FACE_RESTR_INTERIOR)
-
- ambient_dim = mesh.ambient_dim
- dim = mesh.dim
-
- sym_surf_normal = sym.project(dv, df)(
- sym.surface_normal(ambient_dim, dim=dim, dd=dv).as_vector()
- )
- sym_surf_normal = sym_surf_normal / sym.sqrt(sum(sym_surf_normal**2))
-
- sym_face_normal_i = sym.normal(df, ambient_dim, dim=dim - 1)
- sym_face_normal_e = sym.OppositeInteriorFaceSwap(df)(sym_face_normal_i)
-
- if mesh.ambient_dim == 3:
- # NOTE: there's only one face tangent in 3d
- sym_face_tangent = (
- sym.pseudoscalar(ambient_dim, dim - 1, dd=df)
- / sym.area_element(ambient_dim, dim - 1, dd=df)).as_vector()
-
- # }}}
-
- # {{{ checks
-
- def _eval_error(x):
- return bind(discr, sym.norm(np.inf, sym.var("x", dd=df), dd=df))(actx, x=x)
-
- rtol = 1.0e-14
-
- surf_normal = bind(discr, sym_surf_normal)(actx)
-
- face_normal_i = bind(discr, sym_face_normal_i)(actx)
- face_normal_e = bind(discr, sym_face_normal_e)(actx)
-
- # check interpolated surface normal is orthogonal to face normal
- error = _eval_error(surf_normal.dot(face_normal_i))
- logger.info("error[n_dot_i]: %.5e", error)
- assert error < rtol
-
- # check angle between two neighboring elements
- error = _eval_error(face_normal_i.dot(face_normal_e) + 1.0)
- logger.info("error[i_dot_e]: %.5e", error)
- assert error > rtol
-
- # check orthogonality with face tangent
- if ambient_dim == 3:
- face_tangent = bind(discr, sym_face_tangent)(actx)
-
- error = _eval_error(face_tangent.dot(face_normal_i))
- logger.info("error[t_dot_i]: %.5e", error)
- assert error < 5 * rtol
-
- # }}}
-
-# }}}
-
-
-# {{{ diff operator
-
-@pytest.mark.parametrize("dim", [1, 2, 3])
-def test_tri_diff_mat(actx_factory, dim, order=4):
- """Check differentiation matrix along the coordinate axes on a disk
-
- Uses sines as the function to differentiate.
- """
-
- actx = actx_factory()
-
- from pytools.convergence import EOCRecorder
- axis_eoc_recs = [EOCRecorder() for axis in range(dim)]
-
- for n in [4, 8, 16]:
- mesh = mgen.generate_regular_rect_mesh(a=(-0.5,)*dim, b=(0.5,)*dim,
- nelements_per_axis=(n,)*dim, order=4)
-
- discr = DiscretizationCollection(actx, mesh, order=4)
- nabla = sym.nabla(dim)
-
- for axis in range(dim):
- x = sym.nodes(dim)
-
- f = bind(discr, sym.sin(3*x[axis]))(actx)
- df = bind(discr, 3*sym.cos(3*x[axis]))(actx)
-
- sym_op = nabla[axis](sym.var("f"))
- bound_op = bind(discr, sym_op)
- df_num = bound_op(f=f)
-
- linf_error = actx.np.linalg.norm(df_num - df, ord=np.inf)
- axis_eoc_recs[axis].add_data_point(1/n, linf_error)
-
- for axis, eoc_rec in enumerate(axis_eoc_recs):
- logger.info("axis %d\n%s", axis, eoc_rec)
- assert eoc_rec.order_estimate() > order - 0.25
-
-# }}}
-
-
-# {{{ divergence theorem
-
-def test_2d_gauss_theorem(actx_factory):
- """Verify Gauss's theorem explicitly on a mesh"""
-
- pytest.importorskip("meshpy")
-
- from meshpy.geometry import make_circle, GeometryBuilder
- from meshpy.triangle import MeshInfo, build
-
- geob = GeometryBuilder()
- geob.add_geometry(*make_circle(1))
- mesh_info = MeshInfo()
- geob.set(mesh_info)
-
- mesh_info = build(mesh_info)
-
- from meshmode.mesh.io import from_meshpy
- from meshmode.mesh import BTAG_ALL
-
- mesh = from_meshpy(mesh_info, order=1)
-
- actx = actx_factory()
-
- discr = DiscretizationCollection(actx, mesh, order=2)
-
- def f(x):
- return flat_obj_array(
- sym.sin(3*x[0])+sym.cos(3*x[1]),
- sym.sin(2*x[0])+sym.cos(x[1]))
-
- gauss_err = bind(discr,
- sym.integral((
- sym.nabla(2) * f(sym.nodes(2))
- ).sum())
- - # noqa: W504
- sym.integral(
- sym.project("vol", BTAG_ALL)(f(sym.nodes(2)))
- .dot(sym.normal(BTAG_ALL, 2)),
- dd=BTAG_ALL)
- )(actx)
-
- assert abs(gauss_err) < 1e-13
-
-
-@pytest.mark.parametrize("mesh_name", ["2-1-ellipse", "spheroid"])
-def test_surface_divergence_theorem(actx_factory, mesh_name, visualize=False):
- r"""Check the surface divergence theorem.
-
- .. math::
-
- \int_Sigma \phi \nabla_i f_i =
- \int_\Sigma \nabla_i \phi f_i +
- \int_\Sigma \kappa \phi f_i n_i +
- \int_{\partial \Sigma} \phi f_i m_i
-
- where :math:`n_i` is the surface normal and :class:`m_i` is the
- face normal (which should be orthogonal to both the surface normal
- and the face tangent).
- """
- actx = actx_factory()
-
- # {{{ cases
-
- if mesh_name == "2-1-ellipse":
- from mesh_data import EllipseMeshBuilder
- builder = EllipseMeshBuilder(radius=3.1, aspect_ratio=2.0)
- elif mesh_name == "spheroid":
- from mesh_data import SpheroidMeshBuilder
- builder = SpheroidMeshBuilder()
- elif mesh_name == "circle":
- from mesh_data import EllipseMeshBuilder
- builder = EllipseMeshBuilder(radius=1.0, aspect_ratio=1.0)
- elif mesh_name == "starfish":
- from mesh_data import StarfishMeshBuilder
- builder = StarfishMeshBuilder()
- elif mesh_name == "sphere":
- from mesh_data import SphereMeshBuilder
- builder = SphereMeshBuilder(radius=1.0, mesh_order=16)
- else:
- raise ValueError("unknown mesh name: %s" % mesh_name)
-
- # }}}
-
- # {{{ convergene
-
- def f(x):
- return flat_obj_array(
- sym.sin(3*x[1]) + sym.cos(3*x[0]) + 1.0,
- sym.sin(2*x[0]) + sym.cos(x[1]),
- 3.0 * sym.cos(x[0] / 2) + sym.cos(x[1]),
- )[:ambient_dim]
-
- from pytools.convergence import EOCRecorder
- eoc_global = EOCRecorder()
- eoc_local = EOCRecorder()
-
- theta = np.pi / 3.33
- ambient_dim = builder.ambient_dim
- if ambient_dim == 2:
- mesh_rotation = np.array([
- [np.cos(theta), -np.sin(theta)],
- [np.sin(theta), np.cos(theta)],
- ])
- else:
- mesh_rotation = np.array([
- [1.0, 0.0, 0.0],
- [0.0, np.cos(theta), -np.sin(theta)],
- [0.0, np.sin(theta), np.cos(theta)],
- ])
-
- mesh_offset = np.array([0.33, -0.21, 0.0])[:ambient_dim]
-
- for i, resolution in enumerate(builder.resolutions):
- from meshmode.mesh.processing import affine_map
- from meshmode.discretization.connection import FACE_RESTR_ALL
-
- mesh = builder.get_mesh(resolution, builder.mesh_order)
- mesh = affine_map(mesh, A=mesh_rotation, b=mesh_offset)
-
- from meshmode.discretization.poly_element import \
- QuadratureSimplexGroupFactory
- discr = DiscretizationCollection(
- actx, mesh, order=builder.order,
- discr_tag_to_group_factory={
- "product": QuadratureSimplexGroupFactory(2 * builder.order)
- }
- )
-
- volume = discr.discr_from_dd(dof_desc.DD_VOLUME)
- logger.info("ndofs: %d", volume.ndofs)
- logger.info("nelements: %d", volume.mesh.nelements)
-
- dd = dof_desc.DD_VOLUME
- dq = dd.with_discr_tag("product")
- df = dof_desc.as_dofdesc(FACE_RESTR_ALL)
- ambient_dim = discr.ambient_dim
- dim = discr.dim
-
- # variables
- sym_f = f(sym.nodes(ambient_dim, dd=dd))
- sym_f_quad = f(sym.nodes(ambient_dim, dd=dq))
- sym_kappa = sym.summed_curvature(ambient_dim, dim=dim, dd=dq)
- sym_normal = sym.surface_normal(ambient_dim, dim=dim, dd=dq).as_vector()
-
- sym_face_normal = sym.normal(df, ambient_dim, dim=dim - 1)
- sym_face_f = sym.project(dd, df)(sym_f)
-
- # operators
- sym_stiff = sum(
- sym.StiffnessOperator(d)(f) for d, f in enumerate(sym_f)
- )
- sym_stiff_t = sum(
- sym.StiffnessTOperator(d)(f) for d, f in enumerate(sym_f)
- )
- sym_k = sym.MassOperator(dq, dd)(sym_kappa * sym_f_quad.dot(sym_normal))
- sym_flux = sym.FaceMassOperator()(sym_face_f.dot(sym_face_normal))
-
- # sum everything up
- sym_op_global = sym.NodalSum(dd)(
- sym_stiff - (sym_stiff_t + sym_k))
- sym_op_local = sym.ElementwiseSumOperator(dd)(
- sym_stiff - (sym_stiff_t + sym_k + sym_flux))
-
- # evaluate
- op_global = bind(discr, sym_op_global)(actx)
- op_local = bind(discr, sym_op_local)(actx)
-
- err_global = abs(op_global)
- err_local = bind(discr, sym.norm(np.inf, sym.var("x")))(actx, x=op_local)
- logger.info("errors: global %.5e local %.5e", err_global, err_local)
-
- # compute max element size
- h_max = bind(discr, sym.h_max_from_volume(
- discr.ambient_dim, dim=discr.dim, dd=dd))(actx)
- eoc_global.add_data_point(h_max, err_global)
- eoc_local.add_data_point(h_max, err_local)
-
- if visualize:
- from grudge.shortcuts import make_visualizer
- vis = make_visualizer(discr, vis_order=builder.order)
-
- filename = f"surface_divergence_theorem_{mesh_name}_{i:04d}.vtu"
- vis.write_vtk_file(filename, [
- ("r", actx.np.log10(op_local))
- ], overwrite=True)
-
- # }}}
-
- order = min(builder.order, builder.mesh_order) - 0.5
- logger.info("\n%s", str(eoc_global))
- logger.info("\n%s", str(eoc_local))
-
- assert eoc_global.max_error() < 1.0e-12 \
- or eoc_global.order_estimate() > order - 0.5
-
- assert eoc_local.max_error() < 1.0e-12 \
- or eoc_local.order_estimate() > order - 0.5
-
-# }}}
-
-
-# {{{ operator collector determinism
-
-def test_op_collector_order_determinism():
- class TestOperator(sym.Operator):
-
- def __init__(self):
- sym.Operator.__init__(self, dof_desc.DD_VOLUME, dof_desc.DD_VOLUME)
-
- mapper_method = "map_test_operator"
-
- from grudge.symbolic.mappers import BoundOperatorCollector
-
- class TestBoundOperatorCollector(BoundOperatorCollector):
-
- def map_test_operator(self, expr):
- return self.map_operator(expr)
-
- v0 = sym.var("v0")
- ob0 = sym.OperatorBinding(TestOperator(), v0)
-
- v1 = sym.var("v1")
- ob1 = sym.OperatorBinding(TestOperator(), v1)
-
- # The output order isn't significant, but it should always be the same.
- assert list(TestBoundOperatorCollector(TestOperator)(ob0 + ob1)) == [ob0, ob1]
-
-# }}}
-
-
-# {{{ bessel
-
-def test_bessel(actx_factory):
- actx = actx_factory()
-
- dims = 2
-
- mesh = mgen.generate_regular_rect_mesh(
- a=(0.1,)*dims,
- b=(1.0,)*dims,
- nelements_per_axis=(8,)*dims)
-
- discr = DiscretizationCollection(actx, mesh, order=3)
-
- nodes = sym.nodes(dims)
- r = sym.cse(sym.sqrt(nodes[0]**2 + nodes[1]**2))
-
- # https://dlmf.nist.gov/10.6.1
- n = 3
- bessel_zero = (
- sym.bessel_j(n+1, r)
- + sym.bessel_j(n-1, r)
- - 2*n/r * sym.bessel_j(n, r))
-
- z = bind(discr, sym.norm(2, bessel_zero))(actx)
-
- assert z < 1e-15
-
-# }}}
-
-
-# {{{ function symbol
-
-def test_external_call(actx_factory):
- actx = actx_factory()
-
- def double(queue, x):
- return 2 * x
-
- dims = 2
-
- mesh = mgen.generate_regular_rect_mesh(
- a=(0,) * dims, b=(1,) * dims, nelements_per_axis=(4,) * dims)
- discr = DiscretizationCollection(actx, mesh, order=1)
-
- ones = sym.Ones(dof_desc.DD_VOLUME)
- op = (
- ones * 3
- + sym.FunctionSymbol("double")(ones))
-
- from grudge.function_registry import (
- base_function_registry, register_external_function)
-
- freg = register_external_function(
- base_function_registry,
- "double",
- implementation=double,
- dd=dof_desc.DD_VOLUME)
-
- bound_op = bind(discr, op, function_registry=freg)
-
- result = bound_op(actx, double=double)
- assert actx.to_numpy(flatten(result) == 5).all()
-
-
-@pytest.mark.parametrize("array_type", ["scalar", "vector"])
-def test_function_symbol_array(actx_factory, array_type):
- """Test if `FunctionSymbol` distributed properly over object arrays."""
-
- actx = actx_factory()
-
- dim = 2
- mesh = mgen.generate_regular_rect_mesh(
- a=(-0.5,)*dim, b=(0.5,)*dim,
- nelements_per_axis=(8,)*dim, order=4)
- discr = DiscretizationCollection(actx, mesh, order=4)
- volume_discr = discr.discr_from_dd(dof_desc.DD_VOLUME)
-
- if array_type == "scalar":
- sym_x = sym.var("x")
- x = thaw(actx, actx.np.cos(volume_discr.nodes()[0]))
- elif array_type == "vector":
- sym_x = sym.make_sym_array("x", dim)
- x = thaw(actx, volume_discr.nodes())
- else:
- raise ValueError("unknown array type")
-
- norm = bind(discr, sym.norm(2, sym_x))(x=x)
- assert isinstance(norm, float)
-
-# }}}
-
-
-@pytest.mark.parametrize("p", [2, np.inf])
-def test_norm_obj_array(actx_factory, p):
- """Test :func:`grudge.symbolic.operators.norm` for object arrays."""
-
- actx = actx_factory()
-
- dim = 2
- mesh = mgen.generate_regular_rect_mesh(
- a=(-0.5,)*dim, b=(0.5,)*dim,
- nelements_per_axis=(8,)*dim, order=1)
- discr = DiscretizationCollection(actx, mesh, order=4)
-
- w = make_obj_array([1.0, 2.0, 3.0])[:dim]
-
- # {{ scalar
-
- sym_w = sym.var("w")
- norm = bind(discr, sym.norm(p, sym_w))(actx, w=w[0])
-
- norm_exact = w[0]
- logger.info("norm: %.5e %.5e", norm, norm_exact)
- assert abs(norm - norm_exact) < 1.0e-14
-
- # }}}
-
- # {{{ vector
-
- sym_w = sym.make_sym_array("w", dim)
- norm = bind(discr, sym.norm(p, sym_w))(actx, w=w)
-
- norm_exact = np.sqrt(np.sum(w**2)) if p == 2 else np.max(w)
- logger.info("norm: %.5e %.5e", norm, norm_exact)
- assert abs(norm - norm_exact) < 1.0e-14
-
- # }}}
-
-
-def test_map_if(actx_factory):
- """Test :meth:`grudge.symbolic.execution.ExecutionMapper.map_if` handling
- of scalar conditions.
- """
-
- actx = actx_factory()
-
- dim = 2
- mesh = mgen.generate_regular_rect_mesh(
- a=(-0.5,)*dim, b=(0.5,)*dim,
- nelements_per_axis=(8,)*dim, order=4)
- discr = DiscretizationCollection(actx, mesh, order=4)
-
- sym_if = sym.If(sym.Comparison(2.0, "<", 1.0e-14), 1.0, 2.0)
- bind(discr, sym_if)(actx)
-
-
-def test_empty_boundary(actx_factory):
- # https://github.com/inducer/grudge/issues/54
-
- from meshmode.mesh import BTAG_NONE
-
- actx = actx_factory()
-
- dim = 2
- mesh = mgen.generate_regular_rect_mesh(
- a=(-0.5,)*dim, b=(0.5,)*dim,
- nelements_per_axis=(8,)*dim, order=4)
- discr = DiscretizationCollection(actx, mesh, order=4)
- normal = bind(discr,
- sym.normal(BTAG_NONE, dim, dim=dim - 1))(actx)
- from meshmode.dof_array import DOFArray
- for component in normal:
- assert isinstance(component, DOFArray)
- assert len(component) == len(discr.discr_from_dd(BTAG_NONE).groups)
-
-
-def test_operator_compiler_overwrite(actx_factory):
- """Tests that the same expression in ``eval_code`` and ``discr_code``
- does not confuse the OperatorCompiler in grudge/symbolic/compiler.py.
- """
-
- actx = actx_factory()
-
- ambient_dim = 2
- target_order = 4
-
- from meshmode.mesh.generation import generate_regular_rect_mesh
- mesh = generate_regular_rect_mesh(
- a=(-0.5,)*ambient_dim, b=(0.5,)*ambient_dim,
- n=(8,)*ambient_dim, order=1)
- discr = DiscretizationCollection(actx, mesh, order=target_order)
-
- # {{{ test
-
- sym_u = sym.nodes(ambient_dim)
- sym_div_u = sum(d(u) for d, u in zip(sym.nabla(ambient_dim), sym_u))
-
- div_u = bind(discr, sym_div_u)(actx)
- error = bind(discr, sym.norm(2, sym.var("x")))(actx, x=div_u - discr.dim)
- logger.info("error: %.5e", error)
-
- # }}}
-
-
-@pytest.mark.parametrize("ambient_dim", [
- 2,
- # FIXME, cf. https://github.com/inducer/grudge/pull/78/
- pytest.param(3, marks=pytest.mark.xfail)
- ])
-def test_incorrect_assignment_aggregation(actx_factory, ambient_dim):
- """Tests that the greedy assignemnt aggregation code works on a non-trivial
- expression (on which it didn't work at the time of writing).
- """
-
- actx = actx_factory()
-
- target_order = 4
-
- from meshmode.mesh.generation import generate_regular_rect_mesh
- mesh = generate_regular_rect_mesh(
- a=(-0.5,)*ambient_dim, b=(0.5,)*ambient_dim,
- n=(8,)*ambient_dim, order=1)
- discr = DiscretizationCollection(actx, mesh, order=target_order)
-
- # {{{ test with a relative norm
-
- from grudge.dof_desc import DD_VOLUME
- dd = DD_VOLUME
- sym_x = sym.make_sym_array("y", ambient_dim, dd=dd)
- sym_y = sym.make_sym_array("y", ambient_dim, dd=dd)
-
- sym_norm_y = sym.norm(2, sym_y, dd=dd)
- sym_norm_d = sym.norm(2, sym_x - sym_y, dd=dd)
- sym_op = sym_norm_d / sym_norm_y
- logger.info("%s", sym.pretty(sym_op))
-
- # FIXME: this shouldn't raise a RuntimeError
- with pytest.raises(RuntimeError):
- bind(discr, sym_op)(actx, x=1.0, y=discr.discr_from_dd(dd).nodes())
-
- # }}}
-
- # {{{ test with repeated mass inverses
-
- sym_minv_y = sym.cse(sym.InverseMassOperator()(sym_y), "minv_y")
-
- sym_u = make_obj_array([0.5 * sym.Ones(dd), 0.0, 0.0])[:ambient_dim]
- sym_div_u = sum(d(u) for d, u in zip(sym.nabla(ambient_dim), sym_u))
-
- sym_op = sym.MassOperator(dd)(sym_u) \
- + sym.MassOperator(dd)(sym_minv_y * sym_div_u)
- logger.info("%s", sym.pretty(sym_op))
-
- # FIXME: this shouldn't raise a RuntimeError either
- bind(discr, sym_op)(actx, y=discr.discr_from_dd(dd).nodes())
-
- # }}}
-
-
-# You can test individual routines by typing
-# $ python test_grudge.py 'test_routine()'
-
-if __name__ == "__main__":
- import sys
- if len(sys.argv) > 1:
- exec(sys.argv[1])
- else:
- pytest.main([__file__])
-
-# vim: fdm=marker
--
GitLab