diff --git a/grudge/__init__.py b/grudge/__init__.py
index e69de29bb2d1d6434b8b29ae775ad8c2e48c5391..38a84cd2918b14d6805bcf1506b0d1d8d40bc145 100644
--- a/grudge/__init__.py
+++ b/grudge/__init__.py
@@ -0,0 +1,26 @@
+from __future__ import division, absolute_import
+
+__copyright__ = "Copyright (C) 2015 Andreas Kloeckner"
+
+__license__ = """
+Permission is hereby granted, free of charge, to any person obtaining a copy
+of this software and associated documentation files (the "Software"), to deal
+in the Software without restriction, including without limitation the rights
+to use, copy, modify, merge, publish, distribute, sublicense, and/or sell
+copies of the Software, and to permit persons to whom the Software is
+furnished to do so, subject to the following conditions:
+
+The above copyright notice and this permission notice shall be included in
+all copies or substantial portions of the Software.
+
+THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR
+IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY,
+FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL THE
+AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER
+LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM,
+OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN
+THE SOFTWARE.
+"""
+
+import grudge.symbolic.primitives as sym # noqa
+import grudge.symbolic.flux.primitives as sym_flux # noqa
diff --git a/grudge/models/wave.py b/grudge/models/wave.py
index 67767565449e333bd15372362bd772fa548fe76d..a3a9aaf4448127c0bf4889f92adecef821966d9f 100644
--- a/grudge/models/wave.py
+++ b/grudge/models/wave.py
@@ -31,6 +31,7 @@ import numpy as np
from grudge.models import HyperbolicOperator
from grudge.second_order import CentralSecondDerivative
from meshmode.mesh import BTAG_ALL, BTAG_NONE
+from grudge import sym, sym_flux
# {{{ constant-velocity
@@ -79,15 +80,12 @@ class StrongWaveOperator(HyperbolicOperator):
self.flux_type = flux_type
def flux(self):
- from grudge.flux import FluxVectorPlaceholder, make_normal
-
dim = self.dimensions
- w = FluxVectorPlaceholder(1+dim)
+ w = sym_flux.FluxVectorPlaceholder(1+dim)
u = w[0]
v = w[1:]
- normal = make_normal(dim)
+ normal = sym_flux.make_normal(dim)
- from grudge.tools import join_fields
flux_weak = join_fields(
np.dot(v.avg, normal),
u.avg * normal)
diff --git a/grudge/symbolic/__init__.py b/grudge/symbolic/__init__.py
index 900214e18a0fb1c79b5fe0f445163c18f71151ac..acd34e3532cba7516fab14dc1a9d2af46ddd649d 100644
--- a/grudge/symbolic/__init__.py
+++ b/grudge/symbolic/__init__.py
@@ -1,7 +1,6 @@
"""Building blocks and mappers for operator expression trees."""
-from __future__ import division
-from __future__ import absolute_import
+from __future__ import division, absolute_import
__copyright__ = "Copyright (C) 2008 Andreas Kloeckner"
@@ -24,9 +23,3 @@ LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM,
OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN
THE SOFTWARE.
"""
-
-
-from grudge.symbolic.primitives import * # noqa
-from grudge.symbolic.operators import * # noqa
-from grudge.symbolic.mappers import * # noqa
-from grudge.symbolic.tools import * # noqa
diff --git a/grudge/symbolic/flux/__init__.py b/grudge/symbolic/flux/__init__.py
new file mode 100644
index 0000000000000000000000000000000000000000..e69de29bb2d1d6434b8b29ae775ad8c2e48c5391
diff --git a/grudge/symbolic/flux/mappers.py b/grudge/symbolic/flux/mappers.py
new file mode 100644
index 0000000000000000000000000000000000000000..1533274ae9041c5741a57f5d7390aab74673285f
--- /dev/null
+++ b/grudge/symbolic/flux/mappers.py
@@ -0,0 +1,188 @@
+"""Building blocks for flux computation. Flux compilation."""
+
+__copyright__ = "Copyright (C) 2007 Andreas Kloeckner"
+
+__license__ = """
+Permission is hereby granted, free of charge, to any person obtaining a copy
+of this software and associated documentation files (the "Software"), to deal
+in the Software without restriction, including without limitation the rights
+to use, copy, modify, merge, publish, distribute, sublicense, and/or sell
+copies of the Software, and to permit persons to whom the Software is
+furnished to do so, subject to the following conditions:
+
+The above copyright notice and this permission notice shall be included in
+all copies or substantial portions of the Software.
+
+THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR
+IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY,
+FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL THE
+AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER
+LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM,
+OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN
+THE SOFTWARE.
+"""
+
+
+import pymbolic.mapper.collector
+import pymbolic.mapper.flattener
+import pymbolic.mapper.substitutor
+import pymbolic.mapper.constant_folder
+import pymbolic.mapper.flop_counter
+
+
+class FluxIdentityMapperMixin(object):
+ def map_field_component(self, expr):
+ return expr
+
+ def map_normal(self, expr):
+ # a leaf
+ return expr
+
+ map_element_jacobian = map_normal
+ map_face_jacobian = map_normal
+ map_element_order = map_normal
+ map_local_mesh_size = map_normal
+ map_c_function = map_normal
+
+ def map_scalar_parameter(self, expr):
+ return expr
+
+
+class FluxIdentityMapper(
+ pymbolic.mapper.IdentityMapper,
+ FluxIdentityMapperMixin):
+ pass
+
+
+class FluxSubstitutionMapper(pymbolic.mapper.substitutor.SubstitutionMapper,
+ FluxIdentityMapperMixin):
+ def map_field_component(self, expr):
+ result = self.subst_func(expr)
+ if result is not None:
+ return result
+ else:
+ return expr
+
+ map_normal = map_field_component
+
+
+class FluxStringifyMapper(pymbolic.mapper.stringifier.StringifyMapper):
+ def map_field_component(self, expr, enclosing_prec):
+ if expr.is_interior:
+ return "Int[%d]" % expr.index
+ else:
+ return "Ext[%d]" % expr.index
+
+ def map_normal(self, expr, enclosing_prec):
+ return "Normal(%d)" % expr.axis
+
+ def map_element_jacobian(self, expr, enclosing_prec):
+ return "ElJac"
+
+ def map_face_jacobian(self, expr, enclosing_prec):
+ return "FJac"
+
+ def map_element_order(self, expr, enclosing_prec):
+ return "N"
+
+ def map_local_mesh_size(self, expr, enclosing_prec):
+ return "h"
+
+
+class PrettyFluxStringifyMapper(
+ pymbolic.mapper.stringifier.CSESplittingStringifyMapperMixin,
+ FluxStringifyMapper):
+ pass
+
+
+class FluxFlattenMapper(pymbolic.mapper.flattener.FlattenMapper,
+ FluxIdentityMapperMixin):
+ pass
+
+
+class FluxDependencyMapper(pymbolic.mapper.dependency.DependencyMapper):
+ def map_field_component(self, expr):
+ return set([expr])
+
+ def map_normal(self, expr):
+ return set()
+
+ map_element_jacobian = map_normal
+ map_face_jacobian = map_normal
+ map_element_order = map_normal
+ map_local_mesh_size = map_normal
+ map_c_function = map_normal
+
+ def map_scalar_parameter(self, expr):
+ return set([expr])
+
+
+class FluxTermCollector(pymbolic.mapper.collector.TermCollector,
+ FluxIdentityMapperMixin):
+ pass
+
+
+class FluxAllDependencyMapper(FluxDependencyMapper):
+ def map_normal(self, expr):
+ return set([expr])
+
+ def map_element_order(self, expr):
+ return set([expr])
+
+ def map_local_mesh_size(self, expr):
+ return set([expr])
+
+
+class FluxNormalizationMapper(pymbolic.mapper.collector.TermCollector,
+ FluxIdentityMapperMixin):
+ def get_dependencies(self, expr):
+ return FluxAllDependencyMapper()(expr)
+
+ def map_constant_flux(self, expr):
+ if expr.local_c == expr.neighbor_c:
+ return expr.local_c
+ else:
+ return expr
+
+
+class FluxCCFMapper(pymbolic.mapper.constant_folder.CommutativeConstantFoldingMapper,
+ FluxIdentityMapperMixin):
+ def is_constant(self, expr):
+ return not bool(FluxAllDependencyMapper()(expr))
+
+
+class FluxFlipper(FluxIdentityMapper):
+ def map_field_component(self, expr):
+ return expr.__class__(expr.index, not expr.is_interior)
+
+ def map_normal(self, expr):
+ return -expr
+
+ def map_element_jacobian(self, expr):
+ return expr.__class__(not expr.is_interior)
+
+ map_element_order = map_element_jacobian
+
+
+class FluxFlopCounter(pymbolic.mapper.flop_counter.FlopCounter):
+ def map_normal(self, expr):
+ return 0
+
+ map_element_jacobian = map_normal
+ map_face_jacobian = map_normal
+ map_element_order = map_normal
+ map_local_mesh_size = map_normal
+
+ def map_field_component(self, expr):
+ return 0
+
+ def map_function_symbol(self, expr):
+ return 1
+
+ def map_c_function(self, expr):
+ return 1
+
+ def map_scalar_parameter(self, expr):
+ return 0
+
+# vim: foldmethod=marker
diff --git a/grudge/symbolic/flux/primitives b/grudge/symbolic/flux/primitives
new file mode 100644
index 0000000000000000000000000000000000000000..11eb9218492868dbc21495e82c21da256e8356e4
--- /dev/null
+++ b/grudge/symbolic/flux/primitives
@@ -0,0 +1,256 @@
+"""Building blocks for flux computation. Flux compilation."""
+
+__copyright__ = "Copyright (C) 2007 Andreas Kloeckner"
+
+__license__ = """
+Permission is hereby granted, free of charge, to any person obtaining a copy
+of this software and associated documentation files (the "Software"), to deal
+in the Software without restriction, including without limitation the rights
+to use, copy, modify, merge, publish, distribute, sublicense, and/or sell
+copies of the Software, and to permit persons to whom the Software is
+furnished to do so, subject to the following conditions:
+
+The above copyright notice and this permission notice shall be included in
+all copies or substantial portions of the Software.
+
+THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR
+IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY,
+FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL THE
+AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER
+LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM,
+OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN
+THE SOFTWARE.
+"""
+
+
+import numpy
+import pymbolic.primitives
+
+
+class Flux(pymbolic.primitives.AlgebraicLeaf):
+ def stringifier(self):
+ from grudge.symbolic.flux.mappers import FluxStringifyMapper
+ return FluxStringifyMapper
+
+
+class FluxScalarParameter(pymbolic.primitives.Variable):
+ def __init__(self, name, is_complex=False):
+ pymbolic.primitives.Variable.__init__(self, name)
+ self.is_complex = is_complex
+
+ def get_mapper_method(self, mapper):
+ return mapper.map_scalar_parameter
+
+
+class FieldComponent(Flux):
+ def __init__(self, index, is_interior):
+ self.index = index
+ self.is_interior = is_interior
+
+ def is_equal(self, other):
+ return (isinstance(other, FieldComponent)
+ and self.index == other.index
+ and self.is_interior == other.is_interior
+ )
+
+ def __getinitargs__(self):
+ return self.index, self.is_interior
+
+ def get_hash(self):
+ return hash((
+ self.__class__,
+ self.index,
+ self.is_interior))
+
+ def get_mapper_method(self, mapper):
+ return mapper.map_field_component
+
+
+class Normal(Flux):
+ def __init__(self, axis):
+ self.axis = axis
+
+ def __getinitargs__(self):
+ return self.axis,
+
+ def is_equal(self, other):
+ return isinstance(other, Normal) and self.axis == other.axis
+
+ def get_hash(self):
+ return hash((
+ self.__class__,
+ self.axis))
+
+ def get_mapper_method(self, mapper):
+ return mapper.map_normal
+
+
+class _StatelessFlux(Flux):
+ def __getinitargs__(self):
+ return ()
+
+ def is_equal(self, other):
+ return isinstance(other, self.__class__)
+
+ def get_hash(self):
+ return hash(self.__class__)
+
+
+class _SidedFlux(Flux):
+ def __init__(self, is_interior=True):
+ self.is_interior = is_interior
+
+ def __getinitargs__(self):
+ return (self.is_interior,)
+
+ def is_equal(self, other):
+ return (isinstance(other, self.__class__)
+ and self.is_interior == other.is_interior)
+
+ def get_hash(self):
+ return hash((self.__class__, self.is_interior))
+
+
+class FaceJacobian(_StatelessFlux):
+ def get_mapper_method(self, mapper):
+ return mapper.map_face_jacobian
+
+
+class ElementJacobian(_SidedFlux):
+ def get_mapper_method(self, mapper):
+ return mapper.map_element_jacobian
+
+
+class ElementOrder(_SidedFlux):
+ def get_mapper_method(self, mapper):
+ return mapper.map_element_order
+
+
+class LocalMeshSize(_StatelessFlux):
+ def get_mapper_method(self, mapper):
+ return mapper.map_local_mesh_size
+
+
+def make_penalty_term(power=1):
+ from pymbolic.primitives import CommonSubexpression
+ return CommonSubexpression(
+ (ElementOrder()**2/LocalMeshSize())**power,
+ "penalty")
+
+
+PenaltyTerm = make_penalty_term
+
+
+class FluxFunctionSymbol(pymbolic.primitives.FunctionSymbol):
+ pass
+
+
+class Abs(FluxFunctionSymbol):
+ arg_count = 1
+
+
+class Max(FluxFunctionSymbol):
+ arg_count = 2
+
+
+class Min(FluxFunctionSymbol):
+ arg_count = 2
+
+flux_abs = Abs()
+flux_max = Max()
+flux_min = Min()
+
+
+def norm(v):
+ return numpy.dot(v, v)**0.5
+
+
+def make_normal(dimensions):
+ return numpy.array([Normal(i) for i in range(dimensions)], dtype=object)
+
+
+class FluxConstantPlaceholder(object):
+ def __init__(self, constant):
+ self.constant = constant
+
+ @property
+ def int(self):
+ return self.constant
+
+ @property
+ def ext(self):
+ return self.constant
+
+ @property
+ def avg(self):
+ return self.constant
+
+
+class FluxZeroPlaceholder(FluxConstantPlaceholder):
+ def __init__(self):
+ FluxConstantPlaceholder.__init__(self, 0)
+
+
+class FluxScalarPlaceholder(object):
+ def __init__(self, component=0):
+ self.component = component
+
+ def __str__(self):
+ return "FSP(%d)" % self.component
+
+ @property
+ def int(self):
+ return FieldComponent(self.component, True)
+
+ @property
+ def ext(self):
+ return FieldComponent(self.component, False)
+
+ @property
+ def avg(self):
+ return 0.5*(self.int+self.ext)
+
+
+class FluxVectorPlaceholder(object):
+ def __init__(self, components=None, scalars=None):
+ if not (components is not None or scalars is not None):
+ raise ValueError("either components or scalars must be specified")
+ if components is not None and scalars is not None:
+ raise ValueError("only one of components and scalars "
+ "may be specified")
+
+ # make them arrays for the better indexing
+ if components:
+ self.scalars = numpy.array([
+ FluxScalarPlaceholder(i)
+ for i in range(components)])
+ else:
+ self.scalars = numpy.array(scalars)
+
+ def __len__(self):
+ return len(self.scalars)
+
+ def __str__(self):
+ return "%s(%s)" % (self.__class__.__name__, self.scalars)
+
+ def __getitem__(self, idx):
+ if isinstance(idx, int):
+ return self.scalars[idx]
+ else:
+ return FluxVectorPlaceholder(scalars=self.scalars.__getitem__(idx))
+
+ @property
+ def int(self):
+ return numpy.array([scalar.int for scalar in self.scalars])
+
+ @property
+ def ext(self):
+ return numpy.array([scalar.ext for scalar in self.scalars])
+
+ @property
+ def avg(self):
+ return numpy.array([scalar.avg for scalar in self.scalars])
+
+# }}}
+
+# vim: foldmethod=marker
diff --git a/grudge/symbolic/flux/tools.py b/grudge/symbolic/flux/tools.py
new file mode 100644
index 0000000000000000000000000000000000000000..baaee7d5960e1fca2e4cc18e4ec689f42bc98220
--- /dev/null
+++ b/grudge/symbolic/flux/tools.py
@@ -0,0 +1,85 @@
+"""Flux creation helpers."""
+
+from __future__ import division
+
+__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.
+"""
+
+
+
+
+
+def make_lax_friedrichs_flux(wave_speed, state, fluxes, bdry_tags_states_and_fluxes,
+ strong):
+ from pytools.obj_array import join_fields
+ from grudge.flux import make_normal, FluxVectorPlaceholder, flux_max
+
+ n = len(state)
+ d = len(fluxes)
+ normal = make_normal(d)
+ fvph = FluxVectorPlaceholder(len(state)*(1+d)+1)
+
+ wave_speed_ph = fvph[0]
+ state_ph = fvph[1:1+n]
+ fluxes_ph = [fvph[1+i*n:1+(i+1)*n] for i in range(1, d+1)]
+
+ penalty = flux_max(wave_speed_ph.int,wave_speed_ph.ext)*(state_ph.ext-state_ph.int)
+
+ if not strong:
+ num_flux = 0.5*(sum(n_i*(f_i.int+f_i.ext) for n_i, f_i in zip(normal, fluxes_ph))
+ - penalty)
+ else:
+ num_flux = 0.5*(sum(n_i*(f_i.int-f_i.ext) for n_i, f_i in zip(normal, fluxes_ph))
+ + penalty)
+
+ from grudge.optemplate import get_flux_operator
+ flux_op = get_flux_operator(num_flux)
+ int_operand = join_fields(wave_speed, state, *fluxes)
+
+ from grudge.optemplate import BoundaryPair
+ return (flux_op(int_operand)
+ + sum(
+ flux_op(BoundaryPair(int_operand,
+ join_fields(0, bdry_state, *bdry_fluxes), tag))
+ for tag, bdry_state, bdry_fluxes in bdry_tags_states_and_fluxes))
+
+
+
+
+def make_flux_bilinear_form(testee_expr, test_expr):
+ """Create a grudge flux expression for the bilinear form
+
+ .. math::
+
+ \int_{\Gamma} u f(l_i^+,l_i^-),
+
+ where :math:`u` is the (vector or scalar) unknown, given in *testee_epxr*,
+ :math:`f` is a function given by *test_expr* in terms of the Lagrange polynomials
+ :math:`l^_i^{\pm}` against which we are testing. In *test_expr*, :math:`l_i^{\pm}`
+ are represented by :class:`grudge.flux.FieldComponent(-1, True/False)`.
+ :math:`f` is required to be linear in :math:`l_i^{\pm}`.
+ """
+ # keeping this stub to not lose the design work that went into its signature.
+ raise NotImplementedError
+
+
+
diff --git a/grudge/symbolic/mappers/__init__.py b/grudge/symbolic/mappers/__init__.py
new file mode 100644
index 0000000000000000000000000000000000000000..5c88cdbbc92647213c0ab045be82e9ce062f7249
--- /dev/null
+++ b/grudge/symbolic/mappers/__init__.py
@@ -0,0 +1,1662 @@
+"""Operator template mappers."""
+
+from __future__ import division
+
+__copyright__ = "Copyright (C) 2008 Andreas Kloeckner"
+
+__license__ = """
+Permission is hereby granted, free of charge, to any person obtaining a copy
+of this software and associated documentation files (the "Software"), to deal
+in the Software without restriction, including without limitation the rights
+to use, copy, modify, merge, publish, distribute, sublicense, and/or sell
+copies of the Software, and to permit persons to whom the Software is
+furnished to do so, subject to the following conditions:
+
+The above copyright notice and this permission notice shall be included in
+all copies or substantial portions of the Software.
+
+THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR
+IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY,
+FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL THE
+AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER
+LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM,
+OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN
+THE SOFTWARE.
+"""
+
+
+import 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.flop_counter
+from pymbolic.mapper import CSECachingMapperMixin
+
+
+# {{{ mixins
+
+class LocalOpReducerMixin(object):
+ """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(object):
+ """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_nodal_sum = _map_op_base
+ map_nodal_max = _map_op_base
+ map_nodal_min = _map_op_base
+
+ map_diff_base = _map_op_base
+ map_ref_diff_base = _map_op_base
+ map_elementwise_linear = _map_op_base
+ map_flux_base = _map_op_base
+ map_elementwise_max = _map_op_base
+ map_boundarize = _map_op_base
+ map_flux_exchange = _map_op_base
+ map_quad_grid_upsampler = _map_op_base
+ map_quad_int_faces_grid_upsampler = _map_op_base
+ map_quad_bdry_grid_upsampler = _map_op_base
+
+
+class CombineMapperMixin(object):
+ def map_operator_binding(self, expr):
+ return self.combine([self.rec(expr.op), self.rec(expr.field)])
+
+ def map_boundary_pair(self, expr):
+ return self.combine([self.rec(expr.field), self.rec(expr.bfield)])
+
+
+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 expr.__class__(
+ self.rec(expr.op, *args, **kwargs),
+ self.rec(expr.field, *args, **kwargs))
+
+ def map_boundary_pair(self, expr, *args, **kwargs):
+ assert not isinstance(self, BoundOpMapperMixin), \
+ "IdentityMapper instances cannot be combined with " \
+ "the BoundOpMapperMixin"
+
+ return expr.__class__(
+ self.rec(expr.field, *args, **kwargs),
+ self.rec(expr.bfield, *args, **kwargs),
+ expr.tag)
+
+ 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
+
+ def map_scalar_parameter(self, expr, *args, **kwargs):
+ # it's a leaf--no changing children
+ return expr
+
+ map_c_function = map_scalar_parameter
+
+ map_ones = map_scalar_parameter
+ map_node_coordinate_component = map_scalar_parameter
+ map_normal_component = map_elementwise_linear
+ map_jacobian = map_scalar_parameter
+ map_inverse_metric_derivative = map_scalar_parameter
+ map_forward_metric_derivative = map_scalar_parameter
+
+ map_nodal_sum = map_elementwise_linear
+ map_nodal_min = map_elementwise_linear
+ map_nodal_max = map_elementwise_linear
+
+ map_mass_base = map_elementwise_linear
+ map_ref_mass_base = map_elementwise_linear
+ map_diff_base = map_elementwise_linear
+ map_ref_diff_base = map_elementwise_linear
+ map_flux_base = map_elementwise_linear
+ map_elementwise_max = map_elementwise_linear
+ map_boundarize = map_elementwise_linear
+ map_flux_exchange = map_elementwise_linear
+ map_quad_grid_upsampler = map_elementwise_linear
+ map_quad_int_faces_grid_upsampler = map_elementwise_linear
+ map_quad_bdry_grid_upsampler = map_elementwise_linear
+
+
+class BoundOpMapperMixin(object):
+ 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 set([expr])
+ else:
+ return CombineMapperMixin.map_operator_binding(self, expr)
+
+ def map_operator(self, expr):
+ return set()
+
+ def map_scalar_parameter(self, expr):
+ return set([expr])
+
+ def _map_leaf(self, expr):
+ return set()
+
+ map_ones = _map_leaf
+
+ map_node_coordinate_component = _map_leaf
+ map_normal_component = _map_leaf
+
+ map_jacobian = _map_leaf
+ map_forward_metric_derivative = _map_leaf
+ map_inverse_metric_derivative = _map_leaf
+
+
+class FlopCounter(
+ CombineMapperMixin,
+ pymbolic.mapper.flop_counter.FlopCounter):
+ def map_operator_binding(self, expr):
+ return self.rec(expr.field)
+
+ def map_scalar_parameter(self, expr):
+ return 0
+
+ def map_c_function(self, expr):
+ return 1
+
+ def map_ones(self, expr):
+ return 0
+
+ def map_node_coordinate_component(self, expr):
+ return 0
+
+ def map_normal_component(self, expr):
+ return 0
+
+ map_jacobian = map_normal_component
+ map_forward_metric_derivative = map_normal_component
+ map_inverse_metric_derivative = map_normal_component
+
+
+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
+ from grudge.optemplate import Operator, OperatorBinding
+
+ first = expr.children[0]
+ if isinstance(first, 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 OperatorBinding(first, self.rec(prod))
+ else:
+ return first * self.rec(flattened_product(expr.children[1:]))
+
+# }}}
+
+
+# {{{ operator specializer
+
+class OperatorSpecializer(CSECachingMapperMixin, IdentityMapper):
+ """Guided by a typedict obtained through type inference (i.e. by
+ :class:`grudge.optemplate.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.optemplate.mappers.type_inference.TypeInferrer`.
+ """
+ self.typedict = typedict
+
+ map_common_subexpression_uncached = \
+ IdentityMapper.map_common_subexpression
+
+ def map_operator_binding(self, expr):
+ from grudge.optemplate.primitives import BoundaryPair
+
+ from grudge.optemplate.operators import (
+ MassOperator,
+ QuadratureMassOperator,
+ ReferenceMassOperator,
+ ReferenceQuadratureMassOperator,
+
+ StiffnessTOperator,
+ QuadratureStiffnessTOperator,
+ ReferenceStiffnessTOperator,
+ ReferenceQuadratureStiffnessTOperator,
+
+ QuadratureGridUpsampler, QuadratureBoundaryGridUpsampler,
+ FluxOperatorBase, FluxOperator, QuadratureFluxOperator,
+ BoundaryFluxOperator, QuadratureBoundaryFluxOperator,
+ BoundarizeOperator)
+
+ from grudge.optemplate.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, BoundaryPair)
+ has_quad_operand = False
+ else:
+ has_quad_operand = isinstance(field_repr_tag,
+ QuadratureRepresentation)
+ # }}}
+
+ # Based on where this is run in the optemplate processing pipeline,
+ # we may encounter both reference and non-reference operators.
+
+ # {{{ elementwise operators
+ if isinstance(expr.op, MassOperator) and has_quad_operand:
+ return QuadratureMassOperator(
+ field_repr_tag.quadrature_tag)(self.rec(expr.field))
+
+ elif isinstance(expr.op, ReferenceMassOperator) and has_quad_operand:
+ return ReferenceQuadratureMassOperator(
+ field_repr_tag.quadrature_tag)(self.rec(expr.field))
+
+ elif (isinstance(expr.op, StiffnessTOperator) and has_quad_operand):
+ return QuadratureStiffnessTOperator(
+ expr.op.xyz_axis, field_repr_tag.quadrature_tag)(
+ self.rec(expr.field))
+
+ elif (isinstance(expr.op, ReferenceStiffnessTOperator)
+ and has_quad_operand):
+ return ReferenceQuadratureStiffnessTOperator(
+ expr.op.xyz_axis, field_repr_tag.quadrature_tag)(
+ self.rec(expr.field))
+
+ elif (isinstance(expr.op, QuadratureGridUpsampler)
+ and isinstance(field_type, type_info.BoundaryVectorBase)):
+ # potential shortcut:
+ #if (isinstance(expr.field, OperatorBinding)
+ #and isinstance(expr.field.op, BoundarizeOperator)):
+ #return QuadratureBoundarizeOperator(
+ #expr.field.op.tag, expr.op.quadrature_tag)(
+ #self.rec(expr.field.field))
+
+ return QuadratureBoundaryGridUpsampler(
+ expr.op.quadrature_tag, field_type.boundary_tag)(expr.field)
+ # }}}
+
+ elif isinstance(expr.op, BoundarizeOperator) and has_quad_operand:
+ raise TypeError("BoundarizeOperator cannot be applied to "
+ "quadrature-based operands--use QuadUpsample(Boundarize(...))")
+
+ # {{{ flux operator specialization
+ elif isinstance(expr.op, FluxOperatorBase):
+ from pytools.obj_array import with_object_array_or_scalar
+
+ repr_tag_cell = [None]
+
+ def process_flux_arg(flux_arg):
+ arg_repr_tag = self.typedict[flux_arg].repr_tag
+ if repr_tag_cell[0] is None:
+ repr_tag_cell[0] = arg_repr_tag
+ else:
+ # An error for this condition is generated by
+ # the type inference pass.
+
+ assert arg_repr_tag == repr_tag_cell[0]
+
+ is_boundary = isinstance(expr.field, BoundaryPair)
+ if is_boundary:
+ bpair = expr.field
+ with_object_array_or_scalar(process_flux_arg, bpair.field)
+ with_object_array_or_scalar(process_flux_arg, bpair.bfield)
+ else:
+ with_object_array_or_scalar(process_flux_arg, expr.field)
+
+ is_quad = isinstance(repr_tag_cell[0], QuadratureRepresentation)
+ if is_quad:
+ assert not expr.op.is_lift
+ quad_tag = repr_tag_cell[0].quadrature_tag
+
+ new_fld = self.rec(expr.field)
+ flux = expr.op.flux
+
+ if is_boundary:
+ if is_quad:
+ return QuadratureBoundaryFluxOperator(
+ flux, quad_tag, bpair.tag)(new_fld)
+ else:
+ return BoundaryFluxOperator(flux, bpair.tag)(new_fld)
+ else:
+ if is_quad:
+ return QuadratureFluxOperator(flux, quad_tag)(new_fld)
+ else:
+ return FluxOperator(flux, expr.op.is_lift)(new_fld)
+ # }}}
+
+ else:
+ return IdentityMapper.map_operator_binding(self, expr)
+
+ def map_normal_component(self, expr):
+ from grudge.optemplate.mappers.type_inference import (
+ NodalRepresentation)
+
+ expr_type = self.typedict[expr]
+ if not isinstance(
+ expr_type.repr_tag,
+ NodalRepresentation):
+ from grudge.optemplate.primitives import (
+ BoundaryNormalComponent)
+
+ # for now, parts of this are implemented.
+ raise NotImplementedError("normal components on quad. grids")
+
+ return BoundaryNormalComponent(
+ expr.boundary_tag, expr.axis,
+ expr_type.repr_tag.quadrature_tag)
+
+ # a leaf, doesn't change
+ return 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, dimensions):
+ CSECachingMapperMixin.__init__(self)
+ IdentityMapper.__init__(self)
+
+ self.dimensions = dimensions
+
+ map_common_subexpression_uncached = \
+ IdentityMapper.map_common_subexpression
+
+ def map_operator_binding(self, expr):
+ from grudge.optemplate.primitives import (
+ Jacobian, InverseMetricDerivative)
+
+ from grudge.optemplate.operators import (
+ MassOperator,
+ ReferenceMassOperator,
+ QuadratureMassOperator,
+ ReferenceQuadratureMassOperator,
+
+ StiffnessOperator,
+
+ StiffnessTOperator,
+ ReferenceStiffnessTOperator,
+ QuadratureStiffnessTOperator,
+ ReferenceQuadratureStiffnessTOperator,
+
+ InverseMassOperator, ReferenceInverseMassOperator,
+ DifferentiationOperator, ReferenceDifferentiationOperator,
+
+ MInvSTOperator)
+
+ # Global-to-reference is run after operator specialization, so
+ # if we encounter non-quadrature operators here, we know they
+ # must be nodal.
+
+ def rewrite_derivative(ref_class, field, quadrature_tag=None,
+ with_jacobian=True):
+ if quadrature_tag is not None:
+ diff_kwargs = dict(quadrature_tag=quadrature_tag)
+ else:
+ diff_kwargs = {}
+
+ rec_field = self.rec(field)
+ if with_jacobian:
+ rec_field = Jacobian(quadrature_tag) * rec_field
+ return sum(InverseMetricDerivative(None, rst_axis, expr.op.xyz_axis) *
+ ref_class(rst_axis, **diff_kwargs)(rec_field)
+ for rst_axis in range(self.dimensions))
+
+ if isinstance(expr.op, MassOperator):
+ return ReferenceMassOperator()(
+ Jacobian(None) * self.rec(expr.field))
+
+ elif isinstance(expr.op, QuadratureMassOperator):
+ return ReferenceQuadratureMassOperator(
+ expr.op.quadrature_tag)(
+ Jacobian(expr.op.quadrature_tag) * self.rec(expr.field))
+
+ elif isinstance(expr.op, InverseMassOperator):
+ return ReferenceInverseMassOperator()(
+ 1/Jacobian(None) * self.rec(expr.field))
+
+ elif isinstance(expr.op, StiffnessOperator):
+ return ReferenceMassOperator()(Jacobian(None) *
+ self.rec(
+ DifferentiationOperator(expr.op.xyz_axis)(expr.field)))
+
+ elif isinstance(expr.op, DifferentiationOperator):
+ return rewrite_derivative(
+ ReferenceDifferentiationOperator,
+ expr.field, with_jacobian=False)
+
+ elif isinstance(expr.op, StiffnessTOperator):
+ return rewrite_derivative(
+ ReferenceStiffnessTOperator,
+ expr.field)
+
+ elif isinstance(expr.op, QuadratureStiffnessTOperator):
+ return rewrite_derivative(
+ ReferenceQuadratureStiffnessTOperator,
+ expr.field, quadrature_tag=expr.op.quadrature_tag)
+
+ elif isinstance(expr.op, MInvSTOperator):
+ return self.rec(
+ InverseMassOperator()(
+ StiffnessTOperator(expr.op.xyz_axis)(expr.field)))
+
+ else:
+ return IdentityMapper.map_operator_binding(self, expr)
+
+# }}}
+
+
+# {{{ stringification ---------------------------------------------------------
+
+class StringifyMapper(pymbolic.mapper.stringifier.StringifyMapper):
+ def _format_btag(self, tag):
+ from grudge.mesh import SYSTEM_TAGS
+ if tag in SYSTEM_TAGS:
+ return tag.__name__
+ else:
+ return repr(tag)
+
+ def __init__(self, constant_mapper=str, flux_stringify_mapper=None):
+ pymbolic.mapper.stringifier.StringifyMapper.__init__(
+ self, constant_mapper=constant_mapper)
+
+ if flux_stringify_mapper is None:
+ from grudge.flux import FluxStringifyMapper
+ flux_stringify_mapper = FluxStringifyMapper()
+
+ self.flux_stringify_mapper = flux_stringify_mapper
+
+ def map_boundary_pair(self, expr, enclosing_prec):
+ from pymbolic.mapper.stringifier import PREC_NONE
+ return "BPair(%s, %s, %s)" % (
+ self.rec(expr.field, PREC_NONE),
+ self.rec(expr.bfield, PREC_NONE),
+ self._format_btag(expr.tag))
+
+ # {{{ nodal ops
+
+ def map_nodal_sum(self, expr, enclosing_prec):
+ return "NodalSum"
+
+ def map_nodal_max(self, expr, enclosing_prec):
+ return "NodalMax"
+
+ def map_nodal_min(self, expr, enclosing_prec):
+ return "NodalMin"
+
+ # }}}
+
+ # {{{ global differentiation
+
+ def map_diff(self, expr, enclosing_prec):
+ return "Diffx%d" % expr.xyz_axis
+
+ def map_minv_st(self, expr, enclosing_prec):
+ return "MInvSTx%d" % expr.xyz_axis
+
+ def map_stiffness(self, expr, enclosing_prec):
+ return "Stiffx%d" % expr.xyz_axis
+
+ def map_stiffness_t(self, expr, enclosing_prec):
+ return "StiffTx%d" % expr.xyz_axis
+
+ def map_quad_stiffness_t(self, expr, enclosing_prec):
+ return "Q[%s]StiffTx%d" % (
+ expr.quadrature_tag, expr.xyz_axis)
+ # }}}
+
+ # {{{ global mass
+
+ def map_mass(self, expr, enclosing_prec):
+ return "M"
+
+ def map_inverse_mass(self, expr, enclosing_prec):
+ return "InvM"
+
+ def map_quad_mass(self, expr, enclosing_prec):
+ return "Q[%s]M" % expr.quadrature_tag
+
+ # }}}
+
+ # {{{ reference differentiation
+ def map_ref_diff(self, expr, enclosing_prec):
+ return "Diffr%d" % expr.rst_axis
+
+ def map_ref_stiffness_t(self, expr, enclosing_prec):
+ return "StiffTr%d" % expr.rst_axis
+
+ def map_ref_quad_stiffness_t(self, expr, enclosing_prec):
+ return "Q[%s]StiffTr%d" % (
+ expr.quadrature_tag, expr.rst_axis)
+ # }}}
+
+ # {{{ reference mass
+
+ def map_ref_mass(self, expr, enclosing_prec):
+ return "RefM"
+
+ def map_ref_inverse_mass(self, expr, enclosing_prec):
+ return "RefInvM"
+
+ def map_ref_quad_mass(self, expr, enclosing_prec):
+ return "RefQ[%s]M" % expr.quadrature_tag
+
+ # }}}
+
+ def map_elementwise_linear(self, expr, enclosing_prec):
+ return "ElWLin:%s" % expr.__class__.__name__
+
+ # {{{ flux
+
+ def map_flux(self, expr, enclosing_prec):
+ from pymbolic.mapper.stringifier import PREC_NONE
+ return "%s(%s)" % (
+ expr.get_flux_or_lift_text(),
+ self.flux_stringify_mapper(expr.flux, PREC_NONE))
+
+ def map_bdry_flux(self, expr, enclosing_prec):
+ from pymbolic.mapper.stringifier import PREC_NONE
+ return "B[%s]%s(%s)" % (
+ self._format_btag(expr.boundary_tag),
+ expr.get_flux_or_lift_text(),
+ self.flux_stringify_mapper(expr.flux, PREC_NONE))
+
+ def map_quad_flux(self, expr, enclosing_prec):
+ from pymbolic.mapper.stringifier import PREC_NONE
+ return "Q[%s]%s(%s)" % (
+ expr.quadrature_tag,
+ expr.get_flux_or_lift_text(),
+ self.flux_stringify_mapper(expr.flux, PREC_NONE))
+
+ def map_quad_bdry_flux(self, expr, enclosing_prec):
+ from pymbolic.mapper.stringifier import PREC_NONE
+ return "Q[%s]B[%s]%s(%s)" % (
+ expr.quadrature_tag,
+ self._format_btag(expr.boundary_tag),
+ expr.get_flux_or_lift_text(),
+ self.flux_stringify_mapper(expr.flux, PREC_NONE))
+
+ def map_whole_domain_flux(self, expr, enclosing_prec):
+ # used from grudge.backends.cuda.optemplate
+ if expr.is_lift:
+ opname = "WLift"
+ else:
+ opname = "WFlux"
+
+ from pymbolic.mapper.stringifier import PREC_NONE
+ return "%s(%s)" % (opname,
+ self.rec(expr.rebuild_optemplate(), PREC_NONE))
+ # }}}
+
+ def map_elementwise_max(self, expr, enclosing_prec):
+ return "ElWMax"
+
+ def map_boundarize(self, expr, enclosing_prec):
+ return "Boundarize<tag=%s>" % expr.tag
+
+ def map_flux_exchange(self, expr, enclosing_prec):
+ from pymbolic.mapper.stringifier import PREC_NONE
+ return "FExch<idx=%s,rank=%d>(%s)" % (expr.index, expr.rank,
+ ", ".join(self.rec(arg, PREC_NONE) for arg in expr.arg_fields))
+
+ def map_ones(self, expr, enclosing_prec):
+ if expr.quadrature_tag is None:
+ return "Ones"
+ else:
+ return "Q[%s]Ones" % expr.quadrature_tag
+
+ # {{{ geometry data
+
+ def map_node_coordinate_component(self, expr, enclosing_prec):
+ if expr.quadrature_tag is None:
+ return "x[%d]" % expr.axis
+ else:
+ return ("Q[%s]x[%d]" % (expr.quadrature_tag, expr.axis))
+
+ def map_normal_component(self, expr, enclosing_prec):
+ if expr.quadrature_tag is None:
+ return ("Normal<tag=%s>[%d]"
+ % (expr.boundary_tag, expr.axis))
+ else:
+ return ("Q[%s]Normal<tag=%s>[%d]"
+ % (expr.quadrature_tag, expr.boundary_tag, expr.axis))
+
+ def map_jacobian(self, expr, enclosing_prec):
+ if expr.quadrature_tag is None:
+ return "Jac"
+ else:
+ return "JacQ[%s]" % expr.quadrature_tag
+
+ def map_forward_metric_derivative(self, expr, enclosing_prec):
+ result = "dx%d/dr%d" % (expr.xyz_axis, expr.rst_axis)
+ if expr.quadrature_tag is not None:
+ result += "Q[%s]" % expr.quadrature_tag
+ return result
+
+ def map_inverse_metric_derivative(self, expr, enclosing_prec):
+ result = "dr%d/dx%d" % (expr.rst_axis, expr.xyz_axis)
+ if expr.quadrature_tag is not None:
+ result += "Q[%s]" % expr.quadrature_tag
+ return result
+
+ # }}}
+
+ def map_operator_binding(self, expr, enclosing_prec):
+ from pymbolic.mapper.stringifier import PREC_NONE
+ return "<%s>(%s)" % (
+ self.rec(expr.op, PREC_NONE),
+ self.rec(expr.field, PREC_NONE))
+
+ def map_c_function(self, expr, enclosing_prec):
+ return expr.name
+
+ def map_scalar_parameter(self, expr, enclosing_prec):
+ return "ScalarPar[%s]" % expr.name
+
+ def map_quad_grid_upsampler(self, expr, enclosing_prec):
+ return "ToQ[%s]" % expr.quadrature_tag
+
+ def map_quad_int_faces_grid_upsampler(self, expr, enclosing_prec):
+ return "ToIntFaceQ[%s]" % expr.quadrature_tag
+
+ def map_quad_bdry_grid_upsampler(self, expr, enclosing_prec):
+ return "ToBdryQ[%s,%s]" % (expr.quadrature_tag, expr.boundary_tag)
+
+
+class PrettyStringifyMapper(
+ pymbolic.mapper.stringifier.CSESplittingStringifyMapperMixin,
+ StringifyMapper):
+ def __init__(self):
+ pymbolic.mapper.stringifier.CSESplittingStringifyMapperMixin.__init__(self)
+ StringifyMapper.__init__(self)
+
+ self.flux_to_number = {}
+ self.flux_string_list = []
+
+ self.bc_to_number = {}
+ self.bc_string_list = []
+
+ from grudge.flux import PrettyFluxStringifyMapper
+ self.flux_stringify_mapper = PrettyFluxStringifyMapper()
+
+ def get_flux_number(self, flux):
+ try:
+ return self.flux_to_number[flux]
+ except KeyError:
+ from pymbolic.mapper.stringifier import PREC_NONE
+ str_flux = self.flux_stringify_mapper(flux, PREC_NONE)
+
+ flux_number = len(self.flux_to_number)
+ self.flux_string_list.append(str_flux)
+ self.flux_to_number[flux] = flux_number
+ return flux_number
+
+ def map_boundary_pair(self, expr, enclosing_prec):
+ try:
+ bc_number = self.bc_to_number[expr]
+ except KeyError:
+ from pymbolic.mapper.stringifier import PREC_NONE
+ str_bc = StringifyMapper.map_boundary_pair(self, expr, PREC_NONE)
+
+ bc_number = len(self.bc_to_number)
+ self.bc_string_list.append(str_bc)
+ self.bc_to_number[expr] = bc_number
+
+ return "BC%d@%s" % (bc_number, expr.tag)
+
+ def map_operator_binding(self, expr, enclosing_prec):
+ from grudge.optemplate import BoundarizeOperator
+ if isinstance(expr.op, BoundarizeOperator):
+ from pymbolic.mapper.stringifier import PREC_CALL, PREC_SUM
+ return self.parenthesize_if_needed(
+ "%s@%s" % (
+ self.rec(expr.field, PREC_CALL),
+ expr.op.tag),
+ enclosing_prec, PREC_SUM)
+ else:
+ return StringifyMapper.map_operator_binding(
+ self, expr, enclosing_prec)
+
+ def get_bc_strings(self):
+ return ["BC%d : %s" % (i, bc_str)
+ for i, bc_str in enumerate(self.bc_string_list)]
+
+ def get_flux_strings(self):
+ return ["Flux%d : %s" % (i, flux_str)
+ for i, flux_str in enumerate(self.flux_string_list)]
+
+ def map_flux(self, expr, enclosing_prec):
+ return "%s%d" % (
+ expr.get_flux_or_lift_text(),
+ self.get_flux_number(expr.flux))
+
+ def map_bdry_flux(self, expr, enclosing_prec):
+ return "B[%s]%s%d" % (
+ expr.boundary_tag,
+ expr.get_flux_or_lift_text(),
+ self.get_flux_number(expr.flux))
+
+ def map_quad_flux(self, expr, enclosing_prec):
+ return "Q[%s]%s%d" % (
+ expr.quadrature_tag,
+ expr.get_flux_or_lift_text(),
+ self.get_flux_number(expr.flux))
+
+ def map_quad_bdry_flux(self, expr, enclosing_prec):
+ return "Q[%s]B[%s]%s%d" % (
+ expr.quadrature_tag,
+ expr.boundary_tag,
+ expr.get_flux_or_lift_text(),
+ self.get_flux_number(expr.flux))
+
+
+class NoCSEStringifyMapper(StringifyMapper):
+ def map_common_subexpression(self, expr, enclosing_prec):
+ return self.rec(expr.child, enclosing_prec)
+
+# }}}
+
+
+# {{{ quadrature support
+
+class QuadratureUpsamplerRemover(CSECachingMapperMixin, IdentityMapper):
+ def __init__(self, quad_min_degrees, do_warn=True):
+ IdentityMapper.__init__(self)
+ CSECachingMapperMixin.__init__(self)
+ self.quad_min_degrees = quad_min_degrees
+ self.do_warn = do_warn
+
+ map_common_subexpression_uncached = \
+ IdentityMapper.map_common_subexpression
+
+ def map_operator_binding(self, expr):
+ from grudge.optemplate.operators import (
+ QuadratureGridUpsampler,
+ QuadratureInteriorFacesGridUpsampler,
+ QuadratureBoundaryGridUpsampler)
+
+ if isinstance(expr.op, (QuadratureGridUpsampler,
+ QuadratureInteriorFacesGridUpsampler,
+ QuadratureBoundaryGridUpsampler)):
+ try:
+ min_degree = self.quad_min_degrees[expr.op.quadrature_tag]
+ except KeyError:
+ if self.do_warn:
+ from warnings import warn
+ warn("No minimum degree for quadrature tag '%s' specified--"
+ "falling back to nodal evaluation"
+ % expr.op.quadrature_tag)
+ return self.rec(expr.field)
+ else:
+ if min_degree is None:
+ return self.rec(expr.field)
+ else:
+ return expr.op(self.rec(expr.field))
+ else:
+ return IdentityMapper.map_operator_binding(self, expr)
+
+
+class QuadratureDetector(CSECachingMapperMixin, CombineMapper):
+ """For a given expression, this mapper returns the upsampling
+ operator in effect at the root of the expression, or *None*
+ if there isn't one.
+ """
+ class QuadStatusNotKnown:
+ pass
+
+ map_common_subexpression_uncached = \
+ CombineMapper.map_common_subexpression
+
+ def combine(self, values):
+ from pytools import single_valued
+ return single_valued([
+ v for v in values if v is not self.QuadStatusNotKnown])
+
+ def map_variable(self, expr):
+ return None
+
+ def map_constant(self, expr):
+ return self.QuadStatusNotKnown
+
+ def map_operator_binding(self, expr):
+ from grudge.optemplate.operators import (
+ DiffOperatorBase, FluxOperatorBase,
+ MassOperatorBase,
+ QuadratureGridUpsampler,
+ QuadratureInteriorFacesGridUpsampler)
+
+ if isinstance(expr.op, (
+ DiffOperatorBase, FluxOperatorBase,
+ MassOperatorBase)):
+ return None
+ elif isinstance(expr.op, (QuadratureGridUpsampler,
+ QuadratureInteriorFacesGridUpsampler)):
+ return expr.op
+ else:
+ return CombineMapper.map_operator_binding(self, expr)
+
+
+class QuadratureUpsamplerChanger(CSECachingMapperMixin, IdentityMapper):
+ """This mapper descends the expression tree, down to each
+ quadrature-consuming operator (diff, mass) along each branch.
+ It then change
+ """
+ def __init__(self, desired_quad_op):
+ IdentityMapper.__init__(self)
+ CSECachingMapperMixin.__init__(self)
+
+ self.desired_quad_op = desired_quad_op
+
+ map_common_subexpression_uncached = \
+ IdentityMapper.map_common_subexpression
+
+ def map_operator_binding(self, expr):
+ from grudge.optemplate.operators import (
+ DiffOperatorBase, FluxOperatorBase,
+ MassOperatorBase,
+ QuadratureGridUpsampler,
+ QuadratureInteriorFacesGridUpsampler)
+
+ if isinstance(expr.op, (
+ DiffOperatorBase, FluxOperatorBase,
+ MassOperatorBase)):
+ return expr
+ elif isinstance(expr.op, (QuadratureGridUpsampler,
+ QuadratureInteriorFacesGridUpsampler)):
+ return self.desired_quad_op(expr.field)
+ else:
+ return IdentityMapper.map_operator_binding(self, expr)
+
+# }}}
+
+
+# {{{ simplification / optimization
+
+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
+
+ from grudge.optemplate.operators import FluxOperatorBase
+ from grudge.optemplate.primitives import BoundaryPair
+
+ if isinstance(expr.op, FluxOperatorBase):
+ if isinstance(field, BoundaryPair):
+ return self.remove_zeros_from_boundary_flux(expr.op, field)
+ else:
+ return self.remove_zeros_from_interior_flux(expr.op, field)
+ else:
+ return expr.op(field)
+
+ # {{{ remove zeros from interior flux
+ def remove_zeros_from_interior_flux(self, op, vol_field):
+ from pytools.obj_array import is_obj_array
+ if not is_obj_array(vol_field):
+ vol_field = [vol_field]
+
+ from grudge.flux import FieldComponent
+ subst_map = {}
+
+ from grudge.tools import is_zero, make_obj_array
+
+ # process volume field
+ new_vol_field = []
+ new_idx = 0
+ for i, flux_arg in enumerate(vol_field):
+ flux_arg = self.rec(flux_arg)
+
+ if is_zero(flux_arg):
+ subst_map[FieldComponent(i, is_interior=True)] = 0
+ subst_map[FieldComponent(i, is_interior=False)] = 0
+ else:
+ new_vol_field.append(flux_arg)
+ subst_map[FieldComponent(i, is_interior=True)] = \
+ FieldComponent(new_idx, is_interior=True)
+ subst_map[FieldComponent(i, is_interior=False)] = \
+ FieldComponent(new_idx, is_interior=False)
+ new_idx += 1
+
+ # substitute results into flux
+ def sub_flux(expr):
+ return subst_map.get(expr, expr)
+
+ from grudge.flux import FluxSubstitutionMapper
+ new_flux = FluxSubstitutionMapper(sub_flux)(op.flux)
+
+ if is_zero(new_flux):
+ return 0
+ else:
+ return type(op)(new_flux, *op.__getinitargs__()[1:])(
+ make_obj_array(new_vol_field))
+
+ # }}}
+
+ # {{{ remove zeros from boundary flux
+ def remove_zeros_from_boundary_flux(self, op, bpair):
+ vol_field = bpair.field
+ bdry_field = bpair.bfield
+ from pytools.obj_array import is_obj_array
+ if not is_obj_array(vol_field):
+ vol_field = [vol_field]
+ if not is_obj_array(bdry_field):
+ bdry_field = [bdry_field]
+
+ from grudge.flux import FieldComponent
+ subst_map = {}
+
+ # process volume field
+ from grudge.tools import is_zero, make_obj_array
+ new_vol_field = []
+ new_idx = 0
+ for i, flux_arg in enumerate(vol_field):
+ fc = FieldComponent(i, is_interior=True)
+ flux_arg = self.rec(flux_arg)
+
+ if is_zero(flux_arg):
+ subst_map[fc] = 0
+ else:
+ new_vol_field.append(flux_arg)
+ subst_map[fc] = FieldComponent(new_idx, is_interior=True)
+ new_idx += 1
+
+ # process boundary field
+ new_bdry_field = []
+ new_idx = 0
+ for i, flux_arg in enumerate(bdry_field):
+ fc = FieldComponent(i, is_interior=False)
+ flux_arg = self.rec(flux_arg)
+
+ if is_zero(flux_arg):
+ subst_map[fc] = 0
+ else:
+ new_bdry_field.append(flux_arg)
+ subst_map[fc] = FieldComponent(new_idx, is_interior=False)
+ new_idx += 1
+
+ # substitute results into flux
+ def sub_flux(expr):
+ return subst_map.get(expr, expr)
+
+ from grudge.flux import FluxSubstitutionMapper
+ new_flux = FluxSubstitutionMapper(sub_flux)(op.flux)
+
+ if is_zero(new_flux):
+ return 0
+ else:
+ from grudge.optemplate.primitives import BoundaryPair
+ return type(op)(new_flux, *op.__getinitargs__()[1:])(
+ BoundaryPair(
+ make_obj_array(new_vol_field),
+ make_obj_array(new_bdry_field),
+ bpair.tag))
+
+ # }}}
+
+
+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 grudge.optemplate import BoundaryFluxOperatorBase
+
+ if (isinstance(expr.op, BoundaryFluxOperatorBase) and
+ len(self.mesh.tag_to_boundary.get(expr.op.boundary_tag, [])) == 0):
+ return 0
+ else:
+ return IdentityMapper.map_operator_binding(self, expr)
+
+
+class _InnerDerivativeJoiner(pymbolic.mapper.RecursiveMapper):
+ def map_operator_binding(self, expr, derivatives):
+ from grudge.optemplate import DifferentiationOperator
+
+ if isinstance(expr.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.optemplate.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.iteritems():
+ 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_positive(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
+ map_normal_component = map_algebraic_leaf
+ map_jacobian = map_algebraic_leaf
+ map_inverse_metric_derivative = map_algebraic_leaf
+ map_forward_metric_derivative = map_algebraic_leaf
+
+
+class DerivativeJoiner(CSECachingMapperMixin, IdentityMapper):
+ """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.iteritems():
+ derivatives.setdefault(operator, []).extend(operands)
+
+ return result
+
+ derivatives = {}
+ new_children = [invoke_idj(child)
+ for child in expr.children]
+
+ for operator, operands in derivatives.iteritems():
+ 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
+ from grudge.optemplate import InverseMassOperator, OperatorBinding
+
+ if is_zero(expr):
+ return 0
+ else:
+ return OperatorBinding(
+ InverseMassOperator(),
+ self.outer_mass_contractor(expr))
+
+ def map_algebraic_leaf(self, expr):
+ from grudge.optemplate import InverseMassOperator, OperatorBinding
+
+ return OperatorBinding(
+ InverseMassOperator(),
+ self.outer_mass_contractor(expr))
+
+ def map_operator_binding(self, binding):
+ from grudge.optemplate import (
+ MassOperator, StiffnessOperator, StiffnessTOperator,
+ DifferentiationOperator,
+ MInvSTOperator, InverseMassOperator,
+ FluxOperator, BoundaryFluxOperator)
+
+ if isinstance(binding.op, MassOperator):
+ return binding.field
+ elif isinstance(binding.op, StiffnessOperator):
+ return DifferentiationOperator(binding.op.xyz_axis)(
+ self.outer_mass_contractor(binding.field))
+ elif isinstance(binding.op, StiffnessTOperator):
+ return MInvSTOperator(binding.op.xyz_axis)(
+ self.outer_mass_contractor(binding.field))
+ elif isinstance(binding.op, FluxOperator):
+ assert not binding.op.is_lift
+
+ return FluxOperator(binding.op.flux, is_lift=True)(
+ self.outer_mass_contractor(binding.field))
+ elif isinstance(binding.op, BoundaryFluxOperator):
+ assert not binding.op.is_lift
+
+ return BoundaryFluxOperator(binding.op.flux,
+ binding.op.boundary_tag, is_lift=True)(
+ self.outer_mass_contractor(binding.field))
+ else:
+ self.extra_operator_count += 1
+ return InverseMassOperator()(
+ self.outer_mass_contractor(binding))
+
+ def map_sum(self, expr):
+ return expr.__class__(tuple(self.rec(child) for child in expr.children))
+
+ def map_product(self, expr):
+ from grudge.optemplate import InverseMassOperator, ScalarParameter
+
+ def is_scalar(expr):
+ return isinstance(expr, (int, float, complex, ScalarParameter))
+
+ 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 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_boundary_pair(self, bp):
+ from grudge.optemplate import BoundaryPair
+ return BoundaryPair(self.rec(bp.field), self.rec(bp.bfield), bp.tag)
+
+ def map_operator_binding(self, binding):
+ # we only care about bindings of inverse mass operators
+ from grudge.optemplate import InverseMassOperator
+
+ if isinstance(binding.op, InverseMassOperator):
+ 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):
+ from grudge.optemplate import DiffOperatorBase
+
+ if isinstance(expr.op, DiffOperatorBase):
+ if (self.mesh is not None
+ and expr.op.xyz_axis >= self.mesh.dimensions):
+ 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 optemplate components
+
+class CollectorMixin(OperatorReducerMixin, LocalOpReducerMixin, FluxOpReducerMixin):
+ def combine(self, values):
+ from pytools import flatten
+ return set(flatten(values))
+
+ def map_constant(self, bpair):
+ return set()
+
+ map_operator = map_constant
+ map_variable = map_constant
+ map_ones = map_constant
+ map_node_coordinate_component = map_constant
+ map_normal_component = map_constant
+ map_jacobian = map_constant
+ map_forward_metric_derivative = map_constant
+ map_inverse_metric_derivative = map_constant
+ map_scalar_parameter = map_constant
+ map_c_function = map_constant
+
+ def map_whole_domain_flux(self, expr):
+ result = set()
+
+ for ii in expr.interiors:
+ result.update(self.rec(ii.field_expr))
+
+ for bi in expr.boundaries:
+ result.update(self.rec(bi.bpair.field))
+ result.update(self.rec(bi.bpair.bfield))
+
+ return result
+
+
+class FluxCollector(CSECachingMapperMixin, CollectorMixin, CombineMapper):
+ map_common_subexpression_uncached = \
+ CombineMapper.map_common_subexpression
+
+ def map_operator_binding(self, expr):
+ from grudge.optemplate import FluxOperatorBase
+
+ if isinstance(expr.op, FluxOperatorBase):
+ result = set([expr])
+ else:
+ result = set()
+
+ return result | CombineMapper.map_operator_binding(self, expr)
+
+ def map_whole_domain_flux(self, wdflux):
+ result = set([wdflux])
+
+ for intr in wdflux.interiors:
+ result |= self.rec(intr.field_expr)
+ for bdry in wdflux.boundaries:
+ result |= self.rec(bdry.bpair)
+
+ return result
+
+
+class BoundaryTagCollector(CollectorMixin, CombineMapper):
+ def map_boundary_pair(self, bpair):
+ return set([bpair.tag])
+
+
+class GeometricFactorCollector(CollectorMixin, CombineMapper):
+ def map_jacobian(self, expr):
+ return set([expr])
+
+ map_forward_metric_derivative = map_jacobian
+ map_inverse_metric_derivative = map_jacobian
+
+
+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 = set([expr])
+ else:
+ result = set()
+
+ 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 set([expr])
+
+# }}}
+
+
+# {{{ evaluation
+
+class Evaluator(pymbolic.mapper.evaluator.EvaluationMapper):
+ def map_boundary_pair(self, bp):
+ from grudge.optemplate.primitives import BoundaryPair
+ return BoundaryPair(self.rec(bp.field), self.rec(bp.bfield), bp.tag)
+# }}}
+
+
+# {{{ boundary combiner (used by CUDA backend)
+
+class BoundaryCombiner(CSECachingMapperMixin, IdentityMapper):
+ """Combines inner fluxes and boundary fluxes into a
+ single, whole-domain operator of type
+ :class:`grudge.optemplate.operators.WholeDomainFluxOperator`.
+ """
+ def __init__(self, mesh):
+ self.mesh = mesh
+
+ map_common_subexpression_uncached = \
+ IdentityMapper.map_common_subexpression
+
+ def gather_one_wdflux(self, expressions):
+ from grudge.optemplate.primitives import OperatorBinding, BoundaryPair
+ from grudge.optemplate.operators import WholeDomainFluxOperator
+
+ interiors = []
+ boundaries = []
+ is_lift = None
+ # Since None is a valid value of quadrature tags, use
+ # the empty list to symbolize "not known", and add to
+ # list once something is known.
+ quad_tag = []
+
+ rest = []
+
+ for ch in expressions:
+ from grudge.optemplate.operators import FluxOperatorBase
+ if (isinstance(ch, OperatorBinding)
+ and isinstance(ch.op, FluxOperatorBase)):
+ skip = False
+
+ my_is_lift = ch.op.is_lift
+
+ if is_lift is None:
+ is_lift = my_is_lift
+ else:
+ if is_lift != my_is_lift:
+ skip = True
+
+ from grudge.optemplate.operators import \
+ QuadratureFluxOperatorBase
+
+ if isinstance(ch.op, QuadratureFluxOperatorBase):
+ my_quad_tag = ch.op.quadrature_tag
+ else:
+ my_quad_tag = None
+
+ if quad_tag:
+ if quad_tag[0] != my_quad_tag:
+ skip = True
+ else:
+ quad_tag.append(my_quad_tag)
+
+ if skip:
+ rest.append(ch)
+ continue
+
+ if isinstance(ch.field, BoundaryPair):
+ bpair = self.rec(ch.field)
+ if self.mesh.tag_to_boundary.get(bpair.tag, []):
+ boundaries.append(WholeDomainFluxOperator.BoundaryInfo(
+ flux_expr=ch.op.flux,
+ bpair=bpair))
+ else:
+ interiors.append(WholeDomainFluxOperator.InteriorInfo(
+ flux_expr=ch.op.flux,
+ field_expr=self.rec(ch.field)))
+ else:
+ rest.append(ch)
+
+ if interiors or boundaries:
+ wdf = WholeDomainFluxOperator(
+ is_lift=is_lift,
+ interiors=interiors,
+ boundaries=boundaries,
+ quadrature_tag=quad_tag[0])
+ else:
+ wdf = None
+ return wdf, rest
+
+ def map_operator_binding(self, expr):
+ from grudge.optemplate.operators import FluxOperatorBase
+ if isinstance(expr.op, FluxOperatorBase):
+ wdf, rest = self.gather_one_wdflux([expr])
+ assert not rest
+ return wdf
+ else:
+ return IdentityMapper \
+ .map_operator_binding(self, expr)
+
+ def map_sum(self, expr):
+ # FIXME: With flux joining now in the compiler, this is
+ # probably now unnecessary.
+
+ from pymbolic.primitives import flattened_sum
+
+ result = 0
+ expressions = expr.children
+ while True:
+ wdf, expressions = self.gather_one_wdflux(expressions)
+ if wdf is not None:
+ result += wdf
+ else:
+ return result + flattened_sum(self.rec(r_i) for r_i in expressions)
+# }}}
+
+
+# vim: foldmethod=marker
diff --git a/grudge/symbolic/mappers/bc_to_flux.py b/grudge/symbolic/mappers/bc_to_flux.py
new file mode 100644
index 0000000000000000000000000000000000000000..aa2fe090fca46307dad4ff2db2d58634c9dcd927
--- /dev/null
+++ b/grudge/symbolic/mappers/bc_to_flux.py
@@ -0,0 +1,300 @@
+"""Operator template mapper: BC-to-flux rewriting."""
+
+from __future__ import division
+
+__copyright__ = "Copyright (C) 2008 Andreas Kloeckner"
+
+__license__ = """
+Permission is hereby granted, free of charge, to any person obtaining a copy
+of this software and associated documentation files (the "Software"), to deal
+in the Software without restriction, including without limitation the rights
+to use, copy, modify, merge, publish, distribute, sublicense, and/or sell
+copies of the Software, and to permit persons to whom the Software is
+furnished to do so, subject to the following conditions:
+
+The above copyright notice and this permission notice shall be included in
+all copies or substantial portions of the Software.
+
+THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR
+IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY,
+FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL THE
+AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER
+LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM,
+OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN
+THE SOFTWARE.
+"""
+
+
+from pytools import memoize_method
+from pymbolic.mapper import CSECachingMapperMixin
+from grudge.optemplate.mappers import (
+ IdentityMapper, DependencyMapper, CombineMapper,
+ OperatorReducerMixin)
+
+
+class ExpensiveBoundaryOperatorDetector(CombineMapper):
+ def combine(self, values):
+ for val in values:
+ if val:
+ return True
+
+ return False
+
+ def map_operator_binding(self, expr):
+ from grudge.optemplate import (BoundarizeOperator,
+ FluxExchangeOperator,
+ QuadratureGridUpsampler,
+ QuadratureBoundaryGridUpsampler)
+
+ if isinstance(expr.op, BoundarizeOperator):
+ return False
+
+ elif isinstance(expr.op, FluxExchangeOperator):
+ # FIXME: Duplication of these is an even bigger problem!
+ return True
+
+ elif isinstance(expr.op, (
+ QuadratureGridUpsampler,
+ QuadratureBoundaryGridUpsampler)):
+ return True
+
+ else:
+ raise RuntimeError("Found '%s' in a boundary term. "
+ "To the best of my knowledge, no grudge operator applies "
+ "directly to boundary data, so this is likely in error."
+ % expr.op)
+
+ def map_common_subexpression(self, expr):
+ # If there are any expensive operators below here, this
+ # CSE will catch them, so we can easily flux-CSE down to
+ # here.
+
+ return False
+
+ def map_normal_component(self, expr):
+ return False
+
+ map_variable = map_normal_component
+ map_constant = map_normal_component
+
+ @memoize_method
+ def __call__(self, expr):
+ return CombineMapper.__call__(self, expr)
+
+
+class BCToFluxRewriter(CSECachingMapperMixin, IdentityMapper):
+ """Operates on :class:`FluxOperator` instances bound to
+ :class:`BoundaryPair`. If the boundary pair's *bfield* is an expression of
+ what's available in the *field*, we can avoid fetching the data for the
+ explicit boundary condition and just substitute the *bfield* expression
+ into the flux. This mapper does exactly that.
+ """
+
+ map_common_subexpression_uncached = \
+ IdentityMapper.map_common_subexpression
+
+ def __init__(self):
+ self.expensive_bdry_op_detector = \
+ ExpensiveBoundaryOperatorDetector()
+
+ def map_operator_binding(self, expr):
+ from grudge.optemplate.operators import FluxOperatorBase
+ from grudge.optemplate.primitives import BoundaryPair
+ from grudge.flux import FluxSubstitutionMapper, FieldComponent
+
+ if not (isinstance(expr.op, FluxOperatorBase)
+ and isinstance(expr.field, BoundaryPair)):
+ return IdentityMapper.map_operator_binding(self, expr)
+
+ bpair = expr.field
+ vol_field = bpair.field
+ bdry_field = bpair.bfield
+ flux = expr.op.flux
+
+ bdry_dependencies = DependencyMapper(
+ include_calls="descend_args",
+ include_operator_bindings=True)(bdry_field)
+
+ vol_dependencies = DependencyMapper(
+ include_operator_bindings=True)(vol_field)
+
+ vol_bdry_intersection = bdry_dependencies & vol_dependencies
+ if vol_bdry_intersection:
+ raise RuntimeError("Variables are being used as both "
+ "boundary and volume quantities: %s"
+ % ", ".join(str(v) for v in vol_bdry_intersection))
+
+ # Step 1: Find maximal flux-evaluable subexpression of boundary field
+ # in given BoundaryPair.
+
+ class MaxBoundaryFluxEvaluableExpressionFinder(
+ IdentityMapper, OperatorReducerMixin):
+
+ def __init__(self, vol_expr_list, expensive_bdry_op_detector):
+ self.vol_expr_list = vol_expr_list
+ self.vol_expr_to_idx = dict((vol_expr, idx)
+ for idx, vol_expr in enumerate(vol_expr_list))
+
+ self.bdry_expr_list = []
+ self.bdry_expr_to_idx = {}
+
+ self.expensive_bdry_op_detector = expensive_bdry_op_detector
+
+ # {{{ expression registration
+ def register_boundary_expr(self, expr):
+ try:
+ return self.bdry_expr_to_idx[expr]
+ except KeyError:
+ idx = len(self.bdry_expr_to_idx)
+ self.bdry_expr_to_idx[expr] = idx
+ self.bdry_expr_list.append(expr)
+ return idx
+
+ def register_volume_expr(self, expr):
+ try:
+ return self.vol_expr_to_idx[expr]
+ except KeyError:
+ idx = len(self.vol_expr_to_idx)
+ self.vol_expr_to_idx[expr] = idx
+ self.vol_expr_list.append(expr)
+ return idx
+
+ # }}}
+
+ # {{{ map_xxx routines
+
+ @memoize_method
+ def map_common_subexpression(self, expr):
+ # Here we need to decide whether this CSE should be turned into
+ # a flux CSE or not. This is a good idea if the transformed
+ # expression only contains "bare" volume or boundary
+ # expressions. However, as soon as an operator is applied
+ # somewhere in the subexpression, the CSE should not be touched
+ # in order to avoid redundant evaluation of that operator.
+ #
+ # Observe that at the time of this writing (Feb 2010), the only
+ # operators that may occur in boundary expressions are
+ # quadrature-related.
+
+ has_expensive_operators = \
+ self.expensive_bdry_op_detector(expr.child)
+
+ if has_expensive_operators:
+ return FieldComponent(
+ self.register_boundary_expr(expr),
+ is_interior=False)
+ else:
+ return IdentityMapper.map_common_subexpression(self, expr)
+
+ def map_normal(self, expr):
+ raise RuntimeError("Your operator template contains a flux normal. "
+ "You may find this confusing, but you can't do that. "
+ "It turns out that you need to use "
+ "grudge.optemplate.make_normal() for normals in boundary "
+ "terms of operator templates.")
+
+ def map_normal_component(self, expr):
+ if expr.boundary_tag != bpair.tag:
+ raise RuntimeError("BoundaryNormalComponent and BoundaryPair "
+ "do not agree about boundary tag: %s vs %s"
+ % (expr.boundary_tag, bpair.tag))
+
+ from grudge.flux import Normal
+ return Normal(expr.axis)
+
+ def map_variable(self, expr):
+ return FieldComponent(
+ self.register_boundary_expr(expr),
+ is_interior=False)
+
+ map_subscript = map_variable
+
+ def map_operator_binding(self, expr):
+ from grudge.optemplate import (BoundarizeOperator,
+ FluxExchangeOperator,
+ QuadratureGridUpsampler,
+ QuadratureBoundaryGridUpsampler)
+
+ if isinstance(expr.op, BoundarizeOperator):
+ if expr.op.tag != bpair.tag:
+ raise RuntimeError("BoundarizeOperator and BoundaryPair "
+ "do not agree about boundary tag: %s vs %s"
+ % (expr.op.tag, bpair.tag))
+
+ return FieldComponent(
+ self.register_volume_expr(expr.field),
+ is_interior=True)
+
+ elif isinstance(expr.op, FluxExchangeOperator):
+ from grudge.mesh import TAG_RANK_BOUNDARY
+ op_tag = TAG_RANK_BOUNDARY(expr.op.rank)
+ if bpair.tag != op_tag:
+ raise RuntimeError("BoundarizeOperator and "
+ "FluxExchangeOperator do not agree about "
+ "boundary tag: %s vs %s"
+ % (op_tag, bpair.tag))
+ return FieldComponent(
+ self.register_boundary_expr(expr),
+ is_interior=False)
+
+ elif isinstance(expr.op, QuadratureBoundaryGridUpsampler):
+ if bpair.tag != expr.op.boundary_tag:
+ raise RuntimeError("BoundarizeOperator "
+ "and QuadratureBoundaryGridUpsampler "
+ "do not agree about boundary tag: %s vs %s"
+ % (expr.op.boundary_tag, bpair.tag))
+ return FieldComponent(
+ self.register_boundary_expr(expr),
+ is_interior=False)
+
+ elif isinstance(expr.op, QuadratureGridUpsampler):
+ # We're invoked before operator specialization, so we may
+ # see these instead of QuadratureBoundaryGridUpsampler.
+ return FieldComponent(
+ self.register_boundary_expr(expr),
+ is_interior=False)
+
+ else:
+ raise RuntimeError("Found '%s' in a boundary term. "
+ "To the best of my knowledge, no grudge operator applies "
+ "directly to boundary data, so this is likely in error."
+ % expr.op)
+
+ def map_flux_exchange(self, expr):
+ return FieldComponent(
+ self.register_boundary_expr(expr),
+ is_interior=False)
+ # }}}
+
+ from grudge.tools import is_obj_array
+ if not is_obj_array(vol_field):
+ vol_field = [vol_field]
+
+ mbfeef = MaxBoundaryFluxEvaluableExpressionFinder(list(vol_field),
+ self.expensive_bdry_op_detector)
+ #from grudge.optemplate.tools import pretty
+ #print pretty(bdry_field)
+ #raw_input("YO")
+ new_bdry_field = mbfeef(bdry_field)
+
+ # Step II: Substitute the new_bdry_field into the flux.
+ def sub_bdry_into_flux(expr):
+ if isinstance(expr, FieldComponent) and not expr.is_interior:
+ if expr.index == 0 and not is_obj_array(bdry_field):
+ return new_bdry_field
+ else:
+ return new_bdry_field[expr.index]
+ else:
+ return None
+
+ new_flux = FluxSubstitutionMapper(sub_bdry_into_flux)(flux)
+
+ from grudge.tools import is_zero, make_obj_array
+ if is_zero(new_flux):
+ return 0
+ else:
+ return type(expr.op)(new_flux, *expr.op.__getinitargs__()[1:])(
+ BoundaryPair(
+ make_obj_array([self.rec(e) for e in mbfeef.vol_expr_list]),
+ make_obj_array([self.rec(e) for e in mbfeef.bdry_expr_list]),
+ bpair.tag))
diff --git a/grudge/symbolic/mappers/type_inference.py b/grudge/symbolic/mappers/type_inference.py
new file mode 100644
index 0000000000000000000000000000000000000000..a88cfe5a50ec144782189e56f4854525b7e52ebf
--- /dev/null
+++ b/grudge/symbolic/mappers/type_inference.py
@@ -0,0 +1,731 @@
+"""Operator template type inference."""
+
+from __future__ import division
+
+__copyright__ = "Copyright (C) 2008 Andreas Kloeckner"
+
+__license__ = """
+Permission is hereby granted, free of charge, to any person obtaining a copy
+of this software and associated documentation files (the "Software"), to deal
+in the Software without restriction, including without limitation the rights
+to use, copy, modify, merge, publish, distribute, sublicense, and/or sell
+copies of the Software, and to permit persons to whom the Software is
+furnished to do so, subject to the following conditions:
+
+The above copyright notice and this permission notice shall be included in
+all copies or substantial portions of the Software.
+
+THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR
+IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY,
+FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL THE
+AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER
+LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM,
+OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN
+THE SOFTWARE.
+"""
+
+import pymbolic.mapper
+
+
+# {{{ representation tags
+
+class NodalRepresentation(object):
+ """A tag representing nodal representation.
+
+ Volume and boundary vectors below are represented either nodally or on a
+ quadrature grid. This tag expresses one of the two.
+ """
+ def __repr__(self):
+ return "Nodal"
+
+ def __eq__(self, other):
+ return type(self) == type(other)
+
+ def __ne__(self, other):
+ return not self.__eq__(other)
+
+
+class QuadratureRepresentation(object):
+ """A tag representing representation on a quadrature grid tagged with
+ *quadrature_tag".
+
+ Volume and boundary vectors below are represented either nodally or on a
+ quadrature grid. This tag expresses one of the two.
+ """
+ def __init__(self, quadrature_tag):
+ self.quadrature_tag = quadrature_tag
+
+ def __eq__(self, other):
+ return (type(self) == type(other)
+ and self.quadrature_tag == other.quadrature_tag)
+
+ def __ne__(self, other):
+ return not self.__eq__(other)
+
+ def __repr__(self):
+ return "Quadrature(%r)" % self.quadrature_tag
+
+# }}}
+
+
+# {{{ type information --------------------------------------------------------
+
+class type_info:
+ """These classes represent various bits and pieces of information that
+ we may deduce about expressions in our optemplate.
+ """
+
+ # serves only as a namespace, thus lower case
+
+ # {{{ generic type info base classes
+ class TypeInfo(object):
+ def unify(self, other, expr=None):
+ """Return a type that can represent both *self* and *other*.
+ If impossible, raise :exc:`TypeError`. Subtypes should override
+ :meth:`unify_inner`.
+ """
+ # shortcut
+ if self == other:
+ return self
+
+ u_s_o = self.unify_inner(other)
+ u_o_s = other.unify_inner(self)
+
+ if u_s_o is NotImplemented:
+ if u_o_s is NotImplemented:
+ if expr is not None:
+ from grudge.optemplate.tools import pretty
+ raise TypeError("types '%s' and '%s' for '%s' "
+ "cannot be unified" % (self, other,
+ pretty(expr)))
+ else:
+ raise TypeError("types '%s' and '%s' cannot be unified"
+ % (self, other))
+ else:
+ return u_o_s
+ elif u_o_s is NotImplemented:
+ return u_s_o
+
+ if u_s_o != u_o_s:
+ raise RuntimeError("types '%s' and '%s' don't agree about "
+ "their unifier" % (self, other))
+ return u_s_o
+
+ def unify_inner(self, other):
+ """Actual implementation that tries to unify self and other.
+ May return *NotImplemented* to indicate that the reverse unification
+ should be tried. This methods is overriden by derived classes.
+ Derived classes should delegate to base classes if they don't know the
+ answer.
+ """
+ return NotImplemented
+
+ def __eq__(self, other):
+ return (type(self) == type(other)
+ and self.__getinitargs__() == other.__getinitargs__())
+
+ def __ne__(self, other):
+ return not self.__eq__(other)
+
+ class StatelessTypeInfo(TypeInfo):
+ def __getinitargs__(self):
+ return ()
+
+ class FinalType(TypeInfo):
+ """If a :class:`TypeInfo` instance is also an instance of this class,
+ no more information can be added about this type. As a result, this
+ type only unifies with equal instances.
+ """
+ # }}}
+
+ # {{{ simple types: no type, scalar
+ class NoType(StatelessTypeInfo):
+ """Represents "nothing known about type"."""
+ def unify_inner(self, other):
+ return other
+
+ def __repr__(self):
+ return "NoType"
+
+ # this singleton should be the only instance ever created of NoType
+ no_type = NoType()
+
+ class Scalar(StatelessTypeInfo, FinalType):
+ def __repr__(self):
+ return "Scalar"
+ # }}}
+
+ # {{{ tagged type base classes: representation, domain
+ class VectorRepresentationBase(object):
+ def __init__(self, repr_tag):
+ self.repr_tag = repr_tag
+
+ def __getinitargs__(self):
+ return (self.repr_tag,)
+
+ class VolumeVectorBase(object):
+ def __getinitargs__(self):
+ return ()
+
+ class InteriorFacesVectorBase(object):
+ def __getinitargs__(self):
+ return ()
+
+ class BoundaryVectorBase(object):
+ def __init__(self, boundary_tag):
+ self.boundary_tag = boundary_tag
+
+ def __getinitargs__(self):
+ return (self.boundary_tag,)
+ # }}}
+
+ # {{{ single-aspect-known unification helper types
+ class KnownVolume(TypeInfo, VolumeVectorBase):
+ """Type information indicating that this must be a volume vector
+ of unknown representation.
+ """
+
+ def __repr__(self):
+ return "KnownAsVolume"
+
+ def unify_inner(self, other):
+ # Unification with KnownRepresentation is handled in KnownRepresentation.
+ # Here, we only need to unify with VolumeVector.
+
+ if isinstance(other, type_info.VolumeVector):
+ return other
+ else:
+ return type_info.TypeInfo.unify_inner(self, other)
+
+ class KnownInteriorFaces(TypeInfo, InteriorFacesVectorBase):
+ """Type information indicating that this must be a vector
+ of interior face values.
+
+ .. note::
+
+ These vectors are only necessary in a quadrature setting.
+ """
+
+ def __repr__(self):
+ return "KnownAsIntFace"
+
+ def unify_inner(self, other):
+ # Unification with KnownRepresentation is handled in KnownRepresentation.
+ # Here, we only need to unify with InteriorFacesVector.
+
+ if isinstance(other, type_info.InteriorFacesVector):
+ return other
+ elif isinstance(other, type_info.KnownVolume):
+ return type_info.VolumeVector(NodalRepresentation())
+ elif other == type_info.VolumeVector(NodalRepresentation()):
+ return other
+ else:
+ return type_info.TypeInfo.unify_inner(self, other)
+
+ class KnownBoundary(TypeInfo, BoundaryVectorBase):
+ """Type information indicating that this must be a boundary vector."""
+
+ def __repr__(self):
+ return "KnownAsBoundary(%s)" % self.boundary_tag
+
+ def unify_inner(self, other):
+ # Unification with KnownRepresentation is handled in KnownRepresentation.
+ # Here, we only need to unify with BoundaryVector.
+
+ if (isinstance(other, type_info.BoundaryVector)
+ and self.boundary_tag == other.boundary_tag):
+ return other
+ else:
+ return type_info.TypeInfo.unify_inner(self, other)
+
+ class KnownRepresentation(TypeInfo, VectorRepresentationBase):
+ """Type information indicating that the representation (see
+ representation tags, above) is known, but nothing else (e.g. whether
+ this is a boundary or volume vector).
+ """
+ def __repr__(self):
+ return "KnownRepresentation(%s)" % self.repr_tag
+
+ def unify_inner(self, other):
+ if (isinstance(other, type_info.VolumeVector)
+ and self.repr_tag == other.repr_tag):
+ return other
+ elif (isinstance(other, type_info.BoundaryVector)
+ and self.repr_tag == other.repr_tag):
+ return other
+ elif isinstance(other, type_info.KnownVolume):
+ return type_info.VolumeVector(self.repr_tag)
+ elif isinstance(other, type_info.KnownInteriorFaces):
+ if isinstance(self.repr_tag, NodalRepresentation):
+ return type_info.VolumeVector(self.repr_tag)
+ else:
+ return type_info.InteriorFacesVector(self.repr_tag)
+ elif isinstance(other, type_info.KnownBoundary):
+ return type_info.BoundaryVector(other.boundary_tag, self.repr_tag)
+ else:
+ return type_info.TypeInfo.unify_inner(self, other)
+ # }}}
+
+ # {{{ fully specified grudge data types
+ class VolumeVector(FinalType, VectorRepresentationBase, VolumeVectorBase):
+ def __repr__(self):
+ return "Volume(%s)" % self.repr_tag
+
+ class InteriorFacesVector(FinalType, VectorRepresentationBase,
+ InteriorFacesVectorBase):
+ def __init__(self, repr_tag):
+ if not isinstance(repr_tag, QuadratureRepresentation):
+ raise TypeError("InteriorFacesVector is not usable with non-"
+ "quadrature representations")
+ type_info.VectorRepresentationBase.__init__(self, repr_tag)
+
+ def __repr__(self):
+ return "InteriorFaces(%s)" % self.repr_tag
+
+ class BoundaryVector(FinalType, BoundaryVectorBase,
+ VectorRepresentationBase):
+ def __init__(self, boundary_tag, repr_tag):
+ type_info.BoundaryVectorBase.__init__(self, boundary_tag)
+ type_info.VectorRepresentationBase.__init__(self, repr_tag)
+
+ def __repr__(self):
+ return "Boundary(%s, %s)" % (self.boundary_tag, self.repr_tag)
+
+ def __getinitargs__(self):
+ return (self.boundary_tag, self.repr_tag)
+ # }}}
+
+
+# {{{ aspect extraction functions
+
+def extract_representation(ti):
+ try:
+ own_repr_tag = ti.repr_tag
+ except AttributeError:
+ return type_info.no_type
+ else:
+ return type_info.KnownRepresentation(own_repr_tag)
+
+
+def extract_domain(ti):
+ if isinstance(ti, type_info.VolumeVectorBase):
+ return type_info.KnownVolume()
+ elif isinstance(ti, type_info.BoundaryVectorBase):
+ return type_info.KnownBoundary(ti.boundary_tag)
+ else:
+ return type_info.no_type
+
+# }}}
+
+# }}}
+
+
+# {{{ TypeDict helper type
+
+class TypeDict(object):
+ def __init__(self, hints):
+ self.container = hints.copy()
+ self.change_flag = False
+
+ def __getitem__(self, expr):
+ try:
+ return self.container[expr]
+ except KeyError:
+ return type_info.no_type
+
+ def __setitem__(self, expr, new_tp):
+ if new_tp is type_info.no_type:
+ return
+
+ try:
+ old_tp = self.container[expr]
+ except KeyError:
+ self.container[expr] = new_tp
+ self.change_flag = True
+ else:
+ tp = old_tp.unify(new_tp, expr)
+ if tp != old_tp:
+ self.change_flag = True
+ self.container[expr] = tp
+
+ def iteritems(self):
+ return self.container.iteritems()
+
+# }}}
+
+
+# {{{ type inference mapper
+
+class TypeInferrer(pymbolic.mapper.RecursiveMapper):
+ def __init__(self):
+ self.cse_last_results = {}
+
+ def __call__(self, expr, type_hints={}):
+ typedict = TypeDict(type_hints)
+
+ while True:
+ typedict.change_flag = False
+
+ def infer_for_expr(expr):
+ tp = pymbolic.mapper.RecursiveMapper.__call__(self, expr, typedict)
+ typedict[expr] = tp
+
+ # Numpy arrays occur either at the top level or in flux
+ # expressions. This code handles the top level case.
+ from pytools.obj_array import with_object_array_or_scalar
+ with_object_array_or_scalar(infer_for_expr, expr)
+
+ if not typedict.change_flag:
+ # nothing has changed any more, type information has 'converged'
+ break
+
+ # check that type inference completed successfully
+ for expr, tp in typedict.iteritems():
+ if not isinstance(tp, type_info.FinalType):
+ raise RuntimeError("type inference was unable to deduce "
+ "complete type information for '%s' (only '%s')"
+ % (expr, tp))
+
+ return typedict
+
+ def rec(self, expr, typedict):
+ tp = pymbolic.mapper.RecursiveMapper.rec(self, expr, typedict)
+ typedict[expr] = tp
+ return tp
+
+ # Information needs to propagate upward (toward the leaves) *and*
+ # downward (toward the roots) in the expression tree.
+
+ # {{{ base cases
+ def infer_for_children(self, expr, typedict, children):
+ # This routine allows scalars among children and treats them as
+ # not type-changing
+
+ tp = typedict[expr]
+
+ non_scalar_exprs = []
+
+ for child in children:
+ if tp is type_info.no_type:
+ tp = self.rec(child, typedict)
+ if isinstance(tp, type_info.Scalar):
+ tp = type_info.no_type
+ else:
+ non_scalar_exprs.append(child)
+ else:
+ other_tp = self.rec(child, typedict)
+
+ if not isinstance(other_tp, type_info.Scalar):
+ non_scalar_exprs.append(child)
+ tp = tp.unify(other_tp, child)
+
+ for child in non_scalar_exprs:
+ typedict[child] = tp
+
+ if not non_scalar_exprs:
+ tp = type_info.Scalar()
+
+ return tp
+
+ # }}}
+
+ def map_sum(self, expr, typedict):
+ return self.infer_for_children(expr, typedict, expr.children)
+
+ def map_product(self, expr, typedict):
+ return self.infer_for_children(expr, typedict, expr.children)
+
+ def map_quotient(self, expr, typedict):
+ return self.infer_for_children(expr, typedict,
+ children=[expr.numerator, expr.denominator])
+
+ def map_power(self, expr, typedict):
+ return self.infer_for_children(expr, typedict,
+ children=[expr.base, expr.exponent])
+
+ def map_if(self, expr, typedict):
+ return self.infer_for_children(expr, typedict,
+ children=[expr.condition, expr.then, expr.else_])
+
+ def map_comparison(self, expr, typedict):
+ return self.infer_for_children(expr, typedict,
+ children=[expr.left, expr.right])
+
+ def map_if_positive(self, expr, typedict):
+ return self.infer_for_children(expr, typedict,
+ children=[expr.criterion, expr.then, expr.else_])
+
+ def map_call(self, expr, typedict):
+ # assumes functions to be non-type-changing
+ return self.infer_for_children(expr, typedict,
+ children=expr.parameters)
+
+ def map_operator_binding(self, expr, typedict):
+ from grudge.optemplate.operators import (
+ NodalReductionOperator,
+
+ DiffOperatorBase,
+ ReferenceDiffOperatorBase,
+
+ MassOperatorBase,
+ ReferenceMassOperatorBase,
+
+ ElementwiseMaxOperator,
+ BoundarizeOperator, FluxExchangeOperator,
+ FluxOperatorBase,
+ QuadratureGridUpsampler, QuadratureBoundaryGridUpsampler,
+ QuadratureInteriorFacesGridUpsampler,
+
+ MassOperator,
+ ReferenceMassOperator,
+ ReferenceQuadratureMassOperator,
+
+ StiffnessTOperator,
+ ReferenceStiffnessTOperator,
+ ReferenceQuadratureStiffnessTOperator,
+
+ ElementwiseLinearOperator)
+
+ if isinstance(expr.op, NodalReductionOperator):
+ typedict[expr.field] = type_info.KnownVolume()
+ self.rec(expr.field, typedict)
+ return type_info.Scalar()
+
+ elif isinstance(expr.op,
+ (ReferenceQuadratureStiffnessTOperator,
+ ReferenceQuadratureMassOperator)):
+ typedict[expr.field] = type_info.VolumeVector(
+ QuadratureRepresentation(expr.op.quadrature_tag))
+ self.rec(expr.field, typedict)
+ return type_info.VolumeVector(NodalRepresentation())
+
+ elif isinstance(expr.op,
+ (ReferenceStiffnessTOperator, StiffnessTOperator)):
+ # stiffness_T can be specialized for quadrature by OperatorSpecializer
+ typedict[expr.field] = type_info.KnownVolume()
+ self.rec(expr.field, typedict)
+ return type_info.VolumeVector(NodalRepresentation())
+
+ elif isinstance(expr.op,
+ (ReferenceMassOperator, MassOperator)):
+ # mass can be specialized for quadrature by OperatorSpecializer
+ typedict[expr.field] = type_info.KnownVolume()
+ self.rec(expr.field, typedict)
+ return type_info.VolumeVector(NodalRepresentation())
+
+ elif isinstance(expr.op, (
+ DiffOperatorBase,
+ ReferenceDiffOperatorBase,
+ MassOperatorBase,
+ ReferenceMassOperatorBase)):
+ # all other operators are purely nodal
+ typedict[expr.field] = type_info.VolumeVector(NodalRepresentation())
+ self.rec(expr.field, typedict)
+ return type_info.VolumeVector(NodalRepresentation())
+
+ elif isinstance(expr.op, ElementwiseMaxOperator):
+ typedict[expr.field] = typedict[expr].unify(
+ type_info.KnownVolume(), expr.field)
+ return self.rec(expr.field, typedict)
+
+ elif isinstance(expr.op, BoundarizeOperator):
+ # upward propagation: argument has same rep tag as result
+ typedict[expr.field] = type_info.KnownVolume().unify(
+ extract_representation(type_info), expr.field)
+
+ self.rec(expr.field, typedict)
+
+ # downward propagation: result has same rep tag as argument
+ return type_info.KnownBoundary(expr.op.tag).unify(
+ extract_representation(typedict[expr.field]), expr)
+
+ elif isinstance(expr.op, FluxExchangeOperator):
+ raise NotImplementedError
+
+ elif isinstance(expr.op, FluxOperatorBase):
+ from pytools.obj_array import with_object_array_or_scalar
+ from grudge.optemplate.primitives import BoundaryPair
+
+ repr_tag_cell = [type_info.no_type]
+
+ def process_vol_flux_arg(flux_arg):
+ typedict[flux_arg] = type_info.KnownInteriorFaces() \
+ .unify(repr_tag_cell[0], flux_arg)
+ repr_tag_cell[0] = extract_representation(
+ self.rec(flux_arg, typedict))
+
+ if isinstance(expr.field, BoundaryPair):
+ def process_bdry_flux_arg(flux_arg):
+ typedict[flux_arg] = type_info.KnownBoundary(bpair.tag) \
+ .unify(repr_tag_cell[0], flux_arg)
+
+ repr_tag_cell[0] = extract_representation(
+ self.rec(flux_arg, typedict))
+
+ bpair = expr.field
+ with_object_array_or_scalar(process_vol_flux_arg, bpair.field)
+ with_object_array_or_scalar(process_bdry_flux_arg, bpair.bfield)
+ else:
+ with_object_array_or_scalar(process_vol_flux_arg, expr.field)
+
+ return type_info.VolumeVector(NodalRepresentation())
+
+ elif isinstance(expr.op, QuadratureGridUpsampler):
+ typedict[expr.field] = extract_domain(typedict[expr])
+ self.rec(expr.field, typedict)
+ return type_info.KnownRepresentation(
+ QuadratureRepresentation(expr.op.quadrature_tag))\
+ .unify(extract_domain(typedict[expr.field]), expr)
+
+ elif isinstance(expr.op, QuadratureInteriorFacesGridUpsampler):
+ typedict[expr.field] = type_info.VolumeVector(
+ NodalRepresentation())
+ self.rec(expr.field, typedict)
+ return type_info.InteriorFacesVector(
+ QuadratureRepresentation(expr.op.quadrature_tag))
+
+ elif isinstance(expr.op, QuadratureBoundaryGridUpsampler):
+ typedict[expr.field] = type_info.BoundaryVector(
+ expr.op.boundary_tag, NodalRepresentation())
+ self.rec(expr.field, typedict)
+ return type_info.BoundaryVector(
+ expr.op.boundary_tag,
+ QuadratureRepresentation(expr.op.quadrature_tag))
+
+ elif isinstance(expr.op, ElementwiseLinearOperator):
+ typedict[expr.field] = type_info.VolumeVector(NodalRepresentation())
+ self.rec(expr.field, typedict)
+ return type_info.VolumeVector(NodalRepresentation())
+
+ else:
+ raise RuntimeError("TypeInferrer doesn't know how to handle '%s'"
+ % expr.op)
+
+ def map_whole_domain_flux(self, expr, typedict):
+ repr_tag_cell = [type_info.no_type]
+
+ def process_vol_flux_arg(flux_arg):
+ typedict[flux_arg] = type_info.KnownInteriorFaces() \
+ .unify(repr_tag_cell[0], flux_arg)
+ repr_tag_cell[0] = extract_representation(
+ self.rec(flux_arg, typedict))
+
+ def process_bdry_flux_arg(flux_arg):
+ typedict[flux_arg] = type_info.KnownBoundary(bpair.tag) \
+ .unify(repr_tag_cell[0], flux_arg)
+
+ repr_tag_cell[0] = extract_representation(
+ self.rec(flux_arg, typedict))
+
+ from pytools.obj_array import with_object_array_or_scalar
+ for int_flux_info in expr.interiors:
+ with_object_array_or_scalar(process_vol_flux_arg,
+ int_flux_info.field_expr)
+
+ for bdry_flux_info in expr.boundaries:
+ bpair = bdry_flux_info.bpair
+ with_object_array_or_scalar(process_vol_flux_arg, bpair.field)
+ with_object_array_or_scalar(process_bdry_flux_arg, bpair.bfield)
+
+ return type_info.VolumeVector(NodalRepresentation())
+
+ def map_flux_exchange(self, expr, typedict):
+ for arg in expr.arg_fields:
+ typedict[arg] = type_info.VolumeVector(NodalRepresentation())
+
+ from grudge.mesh import TAG_RANK_BOUNDARY
+ return type_info.BoundaryVector(
+ TAG_RANK_BOUNDARY(expr.rank),
+ NodalRepresentation())
+
+ def map_constant(self, expr, typedict):
+ return type_info.Scalar().unify(typedict[expr], expr)
+
+ def map_variable(self, expr, typedict):
+ # user-facing variables are nodal
+ return type_info.KnownRepresentation(NodalRepresentation())\
+ .unify(typedict[expr], expr)
+
+ map_subscript = map_variable
+
+ def map_scalar_parameter(self, expr, typedict):
+ return type_info.Scalar().unify(typedict[expr], expr)
+
+ def map_ones(self, expr, typedict):
+ # FIXME: This is a bit dumb. If the quadrature_tag is None,
+ # we don't know whether the expression was specialized
+ # to 'no quadrature' or if it simply does not know yet
+ # whether it will be on a quadrature grid.
+ if expr.quadrature_tag is not None:
+ return (type_info.VolumeVector(
+ QuadratureRepresentation(expr.quadrature_tag))
+ .unify(typedict[expr], expr))
+ else:
+ return (type_info.VolumeVector(NodalRepresentation())
+ .unify(typedict[expr], expr))
+
+ def map_node_coordinate_component(self, expr, typedict):
+ # FIXME: This is a bit dumb. If the quadrature_tag is None,
+ # we don't know whether the expression was specialized
+ # to 'no quadrature' or if it simply does not know yet
+ # whether it will be on a quadrature grid.
+ if expr.quadrature_tag is not None:
+ return (type_info.VolumeVector(
+ QuadratureRepresentation(expr.quadrature_tag))
+ .unify(typedict[expr], expr))
+ else:
+ return (type_info.VolumeVector(NodalRepresentation())
+ .unify(typedict[expr], expr))
+
+ def map_normal_component(self, expr, typedict):
+ # FIXME: This is a bit dumb. If the quadrature_tag is None,
+ # we don't know whether the expression was specialized
+ # to 'no quadrature' or if it simply does not know yet
+ # whether it will be on a quadrature grid.
+
+ if expr.quadrature_tag is not None:
+ return (type_info.BoundaryVector(expr.boundary_tag,
+ QuadratureRepresentation(expr.quadrature_tag))
+ .unify(typedict[expr], expr))
+ else:
+ return (type_info.KnownBoundary(expr.boundary_tag)
+ .unify(typedict[expr], expr))
+
+ def map_jacobian(self, expr, typedict):
+ return type_info.KnownVolume()
+
+ map_forward_metric_derivative = map_jacobian
+ map_inverse_metric_derivative = map_jacobian
+
+ def map_common_subexpression(self, expr, typedict):
+ outer_tp = typedict[expr]
+
+ last_tp = self.cse_last_results.get(expr, type_info.no_type)
+ if outer_tp != last_tp or last_tp == type_info.no_type:
+ # re-run inner type inference with new outer information
+ typedict[expr.child] = outer_tp
+ new_tp = self.rec(expr.child, typedict)
+
+ # For correct caching, we need to make sure that
+ # information below this level has fully propagated.
+ while True:
+ typedict.change_flag = False
+ new_tp = self.rec(expr.child, typedict)
+
+ if not typedict.change_flag:
+ # nothing has changed any more
+ break
+
+ self.cse_last_results[expr] = new_tp
+ typedict[expr.child] = new_tp
+
+ # we can be sure we *have* changed something
+ typedict.change_flag = True
+ return new_tp
+ else:
+ return last_tp
+
+# }}}
+
+
+# vim: foldmethod=marker
diff --git a/grudge/symbolic/primitives.py b/grudge/symbolic/primitives.py
index e09b707618b93ec4c2d4b4c3e0a84f27ef5ad774..2cfdba75f09e07d92c12f7c8f737fa7ad92521cc 100644
--- a/grudge/symbolic/primitives.py
+++ b/grudge/symbolic/primitives.py
@@ -1,8 +1,6 @@
"""Operator template language: primitives."""
-from __future__ import division
-from __future__ import absolute_import
-from six.moves import range
+from __future__ import division, absolute_import
__copyright__ = "Copyright (C) 2008 Andreas Kloeckner"
@@ -26,10 +24,11 @@ OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN
THE SOFTWARE.
"""
+from six.moves import range
import numpy
import pymbolic.primitives
-import grudge.mesh
+from meshmode.mesh import BTAG_ALL
from pymbolic.primitives import ( # noqa
make_common_subexpression, If, Comparison)