diff --git a/examples/geometry.py b/examples/geometry.py
new file mode 100644
index 0000000000000000000000000000000000000000..f332741576f12174241e6ba706df9a226189fa46
--- /dev/null
+++ b/examples/geometry.py
@@ -0,0 +1,60 @@
+"""Minimal example of a grudge driver."""
+
+from __future__ import division, print_function
+
+__copyright__ = "Copyright (C) 2015 Andreas Kloeckner"
+
+__license__ = """
+Permission is hereby granted, free of charge, to any person obtaining a copy
+of this software and associated documentation files (the "Software"), to deal
+in the Software without restriction, including without limitation the rights
+to use, copy, modify, merge, publish, distribute, sublicense, and/or sell
+copies of the Software, and to permit persons to whom the Software is
+furnished to do so, subject to the following conditions:
+
+The above copyright notice and this permission notice shall be included in
+all copies or substantial portions of the Software.
+
+THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR
+IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY,
+FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL THE
+AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER
+LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM,
+OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN
+THE SOFTWARE.
+"""
+
+import numpy as np # noqa
+import pyopencl as cl
+from grudge import sym, bind, Discretization, shortcuts
+
+
+def main(write_output=True):
+ cl_ctx = cl.create_some_context()
+ queue = cl.CommandQueue(cl_ctx)
+
+ from meshmode.mesh.generation import generate_regular_rect_mesh
+ mesh = generate_regular_rect_mesh(a=(-0.5, -0.5), b=(0.5, 0.5))
+
+ discr = Discretization(cl_ctx, mesh, order=4)
+
+ sym_op = sym.normal(sym.BTAG_ALL, mesh.dim)
+ #sym_op = sym.nodes(mesh.dim, where=sym.BTAG_ALL)
+ print(sym.pretty(sym_op))
+ op = bind(discr, sym_op)
+ print()
+ print(op.code)
+
+ vec = op(queue)
+
+ bvis = shortcuts.make_boundary_visualizer(discr, 4)
+
+ print(vec[0].shape)
+ print(discr.zeros(queue).shape)
+ bvis.write_vtk_file("geo.vtu", [
+ ("normals", vec)
+ ])
+
+
+if __name__ == "__main__":
+ main()
diff --git a/grudge/discretization.py b/grudge/discretization.py
index c9d988448ac8c11043dfa9dc2ad75f3983366b5a..7fb4562df39a35c8c5278f642b47b8fd642bc654 100644
--- a/grudge/discretization.py
+++ b/grudge/discretization.py
@@ -89,5 +89,9 @@ class Discretization(object):
def zeros(self, queue, dtype=None, extra_dims=None):
return self.volume_discr.zeros(queue, dtype, extra_dims=None)
+ def is_volume_where(self, where):
+ from grudge import sym
+ return where == sym.VTAG_ALL
+
# vim: foldmethod=marker
diff --git a/grudge/execution.py b/grudge/execution.py
index df8f0724e116fdddef294b68d764f523a3e787fd..030ab447bd3918347b3df1fde2d835520565015a 100644
--- a/grudge/execution.py
+++ b/grudge/execution.py
@@ -22,8 +22,10 @@ OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN
THE SOFTWARE.
"""
+import six
import numpy as np
import grudge.symbolic.mappers as mappers
+from grudge import sym
import logging
logger = logging.getLogger(__name__)
@@ -56,13 +58,12 @@ class ExecutionMapper(mappers.Evaluator,
if expr.quadrature_tag is not None:
raise NotImplementedError("node coordinate components on quad. grids")
- return self.discr.volume_discr.nodes()[expr.axis] \
- .with_queue(self.queue)
+ if self.discr.is_volume_where(expr.where):
+ discr = self.discr.volume_discr
+ else:
+ discr = self.discr.boundary_discr
- def map_normal_component(self, expr):
- if expr.quadrature_tag is not None:
- raise NotImplementedError("normal components on quad. grids")
- return self.discr.boundary_normals(expr.boundary_tag)[expr.axis]
+ return discr.nodes()[expr.axis].with_queue(self.queue)
def map_boundarize(self, op, field_expr):
return self.discr.boundarize_volume_field(
@@ -72,26 +73,13 @@ class ExecutionMapper(mappers.Evaluator,
def map_scalar_parameter(self, expr):
return self.context[expr.name]
- def map_jacobian(self, expr):
- return self.discr.volume_jacobians(expr.quadrature_tag)
-
- def map_forward_metric_derivative(self, expr):
- return (self.discr.forward_metric_derivatives(expr.quadrature_tag)
- [expr.xyz_axis][expr.rst_axis])
-
- def map_inverse_metric_derivative(self, expr):
- return (self.discr.inverse_metric_derivatives(expr.quadrature_tag)
- [expr.xyz_axis][expr.rst_axis])
-
def map_call(self, expr):
from pymbolic.primitives import Variable
assert isinstance(expr.function, Variable)
- func_name = expr.function.name
- try:
- func = self.discr.exec_functions[func_name]
- except KeyError:
- func = getattr(np, expr.function.name)
+ # FIXME: Make a way to register functions
+ import pyopencl.clmath as clmath
+ func = getattr(clmath, expr.function.name)
return func(*[self.rec(p) for p in expr.parameters])
@@ -241,7 +229,8 @@ class ExecutionMapper(mappers.Evaluator,
# }}}
- # {{{ code execution functions --------------------------------------------
+ # {{{ code execution functions
+
def exec_assign(self, insn):
return [(name, self.rec(expr))
for name, expr in zip(insn.names, insn.exprs)], []
@@ -254,7 +243,8 @@ class ExecutionMapper(mappers.Evaluator,
else:
stats_callback = None
- if insn.flop_count() == 0:
+ # FIXME: Reenable compiled vector exprs
+ if True: # insn.flop_count() == 0:
return [(name, self(expr))
for name, expr in zip(insn.names, insn.exprs)], []
else:
@@ -387,9 +377,22 @@ class ExecutionMapper(mappers.Evaluator,
return result, []
def exec_diff_batch_assign(self, insn):
- rst_diff = self.executor.diff(insn.operators, self.rec(insn.field))
+ field = self.rec(insn.field)
+ repr_op = insn.operators[0]
+ if not isinstance(repr_op, sym.ReferenceDifferentiationOperator):
+ # FIXME
+ raise NotImplementedError()
+
+ if self.discr.is_volume_where(repr_op.where):
+ discr = self.discr.volume_discr
+ else:
+ discr = self.discr.boundary_discr
- return [(name, diff) for name, diff in zip(insn.names, rst_diff)], []
+ return [
+ (name, discr.num_reference_derivative(
+ self.queue, (op.rst_axis,), field)
+ .with_queue(self.queue))
+ for name, op in zip(insn.names, insn.operators)], []
exec_quad_diff_batch_assign = exec_diff_batch_assign
@@ -496,7 +499,21 @@ class Executor(object):
coeffs, matrix, field, out)
def __call__(self, queue, **context):
- return self.code.execute(ExecutionMapper(queue, context, self))
+ import pyopencl.array as cl_array
+
+ def replace_queue(a):
+ if isinstance(a, cl_array.Array):
+ return a.with_queue(queue)
+ else:
+ return a
+
+ from pytools.obj_array import with_object_array_or_scalar
+
+ new_context = {}
+ for name, var in six.iteritems(context):
+ new_context[name] = with_object_array_or_scalar(replace_queue, var)
+
+ return self.code.execute(ExecutionMapper(queue, new_context, self))
# }}}
diff --git a/grudge/shortcuts.py b/grudge/shortcuts.py
index 0de88a50df5527cd9834401f548cbe9f24c2c1c1..0e6940ff6f120018a8b5c3c9af29de04eb1bddcb 100644
--- a/grudge/shortcuts.py
+++ b/grudge/shortcuts.py
@@ -25,6 +25,9 @@ THE SOFTWARE.
"""
+import pyopencl as cl
+
+
def set_up_rk4(field_var_name, dt, fields, rhs, t_start=0):
from leap.rk import LSRK4Method
from dagrt.codegen import PythonCodeGenerator
@@ -39,3 +42,15 @@ def set_up_rk4(field_var_name, dt, fields, rhs, t_start=0):
context={dt_method.component_id: fields})
return dt_stepper
+
+
+def make_visualizer(discr, vis_order):
+ from meshmode.discretization.visualization import make_visualizer
+ with cl.CommandQueue(discr.cl_context) as queue:
+ return make_visualizer(queue, discr.volume_discr, vis_order)
+
+
+def make_boundary_visualizer(discr, vis_order):
+ from meshmode.discretization.visualization import make_visualizer
+ with cl.CommandQueue(discr.cl_context) as queue:
+ return make_visualizer(queue, discr.boundary_discr, vis_order)
diff --git a/grudge/symbolic/compiler.py b/grudge/symbolic/compiler.py
index 2b367255e01eae3236f0f5853502971023646b40..cfbec3c3d7e41e7045bf9c8d157a8312e5821b27 100644
--- a/grudge/symbolic/compiler.py
+++ b/grudge/symbolic/compiler.py
@@ -526,7 +526,7 @@ class Code(object):
if self.static_schedule_attempts:
self.last_schedule = schedule
- from grudge.tools import with_object_array_or_scalar
+ from pytools.obj_array import with_object_array_or_scalar
return with_object_array_or_scalar(exec_mapper, self.result)
# }}}
@@ -642,7 +642,7 @@ class OperatorCompiler(mappers.IdentityMapper):
# {{{ top-level driver ----------------------------------------------------
def __call__(self, expr, type_hints={}):
# Put the result expressions into variables as well.
- expr = sym.make_common_subexpression(expr, "_result")
+ expr = sym.cse(expr, "_result")
from grudge.symbolic.mappers.type_inference import TypeInferrer
self.typedict = TypeInferrer()(expr, type_hints)
@@ -702,7 +702,9 @@ class OperatorCompiler(mappers.IdentityMapper):
# Finally, walk the expression and build the code.
result = super(OperatorCompiler, self).__call__(expr)
- return Code(self.aggregate_assignments(self.code, result), result)
+ # FIXME: Reenable assignment aggregation
+ #return Code(self.aggregate_assignments(self.code, result), result)
+ return Code(self.code, result)
# }}}
diff --git a/grudge/symbolic/mappers/__init__.py b/grudge/symbolic/mappers/__init__.py
index ed8ff2341fed41a0476003eae49f5a5c944e09e2..14fea40a18c1edf788ee660d80a479a12a3e4bfa 100644
--- a/grudge/symbolic/mappers/__init__.py
+++ b/grudge/symbolic/mappers/__init__.py
@@ -186,10 +186,6 @@ class IdentityMapperMixin(LocalOpReducerMixin, FluxOpReducerMixin):
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
@@ -258,11 +254,6 @@ class DependencyMapper(
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(
@@ -283,13 +274,6 @@ class FlopCounter(
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,
@@ -475,27 +459,6 @@ class OperatorSpecializer(CSECachingMapperMixin, IdentityMapper):
else:
return IdentityMapper.map_operator_binding(self, expr)
- def map_normal_component(self, expr):
- from grudge.symbolic.mappers.type_inference import (
- NodalRepresentation)
-
- expr_type = self.typedict[expr]
- if not isinstance(
- expr_type.repr_tag,
- NodalRepresentation):
- from grudge.symbolic.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
-
# }}}
@@ -506,25 +469,32 @@ class GlobalToReferenceMapper(CSECachingMapperMixin, IdentityMapper):
reference elements, together with explicit multiplication by geometric factors.
"""
- def __init__(self, dimensions):
+ def __init__(self, ambient_dim, dim=None):
CSECachingMapperMixin.__init__(self)
IdentityMapper.__init__(self)
- self.dimensions = dimensions
+ if dim is None:
+ dim = ambient_dim
+
+ self.ambient_dim = ambient_dim
+ self.dim = dim
map_common_subexpression_uncached = \
IdentityMapper.map_common_subexpression
def map_operator_binding(self, expr):
- from grudge.symbolic.primitives import (
- Jacobian, InverseMetricDerivative)
-
# 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,
+ jac_notag = sym.area_element(self.ambient_dim, self.dim,
+ where=expr.op.where, quadrature_tag=None)
+
+ def rewrite_derivative(ref_class, field, quadrature_tag, where,
with_jacobian=True):
+ jac_tag = sym.area_element(self.ambient_dim, self.dim,
+ where=expr.op.where, quadrature_tag=quadrature_tag)
+
if quadrature_tag is not None:
diff_kwargs = dict(quadrature_tag=quadrature_tag)
else:
@@ -532,43 +502,50 @@ class GlobalToReferenceMapper(CSECachingMapperMixin, IdentityMapper):
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) *
+ rec_field = jac_tag * rec_field
+ return sum(
+ sym.inverse_metric_derivative(
+ rst_axis, expr.op.xyz_axis,
+ ambient_dim=self.ambient_dim, dim=self.dim) *
ref_class(rst_axis, **diff_kwargs)(rec_field)
for rst_axis in range(self.dimensions))
if isinstance(expr.op, op.MassOperator):
return op.ReferenceMassOperator()(
- Jacobian(None) * self.rec(expr.field))
+ jac_notag * self.rec(expr.field))
elif isinstance(expr.op, op.QuadratureMassOperator):
+ jac_tag = sym.area_element(self.ambient_dim, self.dim,
+ where=expr.op.where,
+ quadrature_tag=expr.op.quadrature_tag)
return op.ReferenceQuadratureMassOperator(
- expr.op.quadrature_tag)(
- Jacobian(expr.op.quadrature_tag) * self.rec(expr.field))
+ expr.op.quadrature_tag, expr.op.where)(
+ jac_tag * self.rec(expr.field))
elif isinstance(expr.op, op.InverseMassOperator):
- return op.ReferenceInverseMassOperator()(
- 1/Jacobian(None) * self.rec(expr.field))
+ return op.ReferenceInverseMassOperator(expr.op.where)(
+ 1/jac_notag * self.rec(expr.field))
elif isinstance(expr.op, op.StiffnessOperator):
- return op.ReferenceMassOperator()(Jacobian(None) *
+ return op.ReferenceMassOperator(expr.op.where)(jac_notag *
self.rec(
op.DifferentiationOperator(expr.op.xyz_axis)(expr.field)))
elif isinstance(expr.op, op.DifferentiationOperator):
return rewrite_derivative(
op.ReferenceDifferentiationOperator,
- expr.field, with_jacobian=False)
+ expr.field, where=expr.op.where, with_jacobian=False)
elif isinstance(expr.op, op.StiffnessTOperator):
return rewrite_derivative(
op.ReferenceStiffnessTOperator,
- expr.field)
+ expr.field, where=expr.op.where)
elif isinstance(expr.op, op.QuadratureStiffnessTOperator):
return rewrite_derivative(
op.ReferenceQuadratureStiffnessTOperator,
- expr.field, quadrature_tag=expr.op.quadrature_tag)
+ expr.field, quadrature_tag=expr.op.quadrature_tag,
+ where=expr.op.where)
elif isinstance(expr.op, op.MInvSTOperator):
return self.rec(
@@ -590,6 +567,18 @@ class StringifyMapper(pymbolic.mapper.stringifier.StringifyMapper):
else:
return repr(tag)
+ def _format_op_props(self, op):
+ result = ""
+
+ if hasattr(op, "where") and op.where is not None:
+ result += "@"+self._format_btag(op)
+
+ if (hasattr(op, "quadrature_tag")
+ and op.quadrature_tag is not None):
+ result += "Q[%s]" % self.quadrature_tag
+
+ return result
+
def __init__(self, constant_mapper=str, flux_stringify_mapper=None):
pymbolic.mapper.stringifier.StringifyMapper.__init__(
self, constant_mapper=constant_mapper)
@@ -623,20 +612,21 @@ class StringifyMapper(pymbolic.mapper.stringifier.StringifyMapper):
# {{{ global differentiation
def map_diff(self, expr, enclosing_prec):
- return "Diffx%d" % expr.xyz_axis
+ return "Diffx%d%s" % (expr.xyz_axi, self._format_op_props(expr))
def map_minv_st(self, expr, enclosing_prec):
- return "MInvSTx%d" % expr.xyz_axis
+ return "MInvSTx%d%s" % (expr.xyz_axi, self._format_op_props(expr))
def map_stiffness(self, expr, enclosing_prec):
- return "Stiffx%d" % expr.xyz_axis
+ return "Stiffx%d%s" % (expr.xyz_axi, self._format_op_props(expr))
def map_stiffness_t(self, expr, enclosing_prec):
- return "StiffTx%d" % expr.xyz_axis
+ return "StiffTx%d%s" % (expr.xyz_axis, self._format_op_props(expr))
def map_quad_stiffness_t(self, expr, enclosing_prec):
- return "Q[%s]StiffTx%d" % (
- expr.quadrature_tag, expr.xyz_axis)
+ return "StiffTx%d%s" % (
+ expr.input_quadrature_tag, expr.xyz_axis,
+ self._format_op_props(expr))
# }}}
# {{{ global mass
@@ -654,31 +644,34 @@ class StringifyMapper(pymbolic.mapper.stringifier.StringifyMapper):
# {{{ reference differentiation
def map_ref_diff(self, expr, enclosing_prec):
- return "Diffr%d" % expr.rst_axis
+ return "Diffr%d%s" % (expr.rst_axis, self._format_op_props(expr))
def map_ref_stiffness_t(self, expr, enclosing_prec):
- return "StiffTr%d" % expr.rst_axis
+ return "StiffTr%d%s" % (expr.rst_axis, self._format_op_props(expr))
def map_ref_quad_stiffness_t(self, expr, enclosing_prec):
- return "Q[%s]StiffTr%d" % (
- expr.quadrature_tag, expr.rst_axis)
+ return "StiffTr%dQ[%s]%s" % (
+ expr.input_quadrature_tag, expr.rst_axis,
+ self._format_op_props(expr))
# }}}
# {{{ reference mass
def map_ref_mass(self, expr, enclosing_prec):
- return "RefM"
+ return "RefM" + self._format_op_props(expr)
def map_ref_inverse_mass(self, expr, enclosing_prec):
- return "RefInvM"
+ return "RefInvM" + self._format_op_props(expr)
def map_ref_quad_mass(self, expr, enclosing_prec):
- return "RefQ[%s]M" % expr.quadrature_tag
+ return "RefM" + self._format_op_props(expr)
# }}}
def map_elementwise_linear(self, expr, enclosing_prec):
- return "ElWLin:%s" % expr.__class__.__name__
+ return "ElWLin:%s%s" % (
+ expr.__class__.__name__,
+ self._format_op_props(expr))
# {{{ flux
@@ -723,10 +716,12 @@ class StringifyMapper(pymbolic.mapper.stringifier.StringifyMapper):
# }}}
def map_elementwise_max(self, expr, enclosing_prec):
- return "ElWMax"
+ return "ElWMax" + self._format_op_props(expr)
def map_boundarize(self, expr, enclosing_prec):
- return "Boundarize<tag=%s>" % self._format_btag(expr.tag)
+ return "Boundarize<tag=%s>%s" % (
+ self._format_btag(expr.tag),
+ self._format_op_props(expr))
def map_flux_exchange(self, expr, enclosing_prec):
from pymbolic.mapper.stringifier import PREC_NONE
@@ -734,10 +729,7 @@ class StringifyMapper(pymbolic.mapper.stringifier.StringifyMapper):
", ".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
+ return "Ones" + self._format_op_props(expr)
# {{{ geometry data
@@ -747,33 +739,6 @@ class StringifyMapper(pymbolic.mapper.stringifier.StringifyMapper):
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]"
- % (self._format_btag(expr.boundary_tag), expr.axis))
- else:
- return ("Q[%s]Normal<tag=%s>[%d]"
- % (expr.quadrature_tag, self._format_btag(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):
@@ -789,14 +754,11 @@ class StringifyMapper(pymbolic.mapper.stringifier.StringifyMapper):
return "ScalarPar[%s]" % expr.name
def map_quad_grid_upsampler(self, expr, enclosing_prec):
- return "ToQ[%s]" % expr.quadrature_tag
+ return "Upsamp" + self._format_op_props(expr)
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,
@@ -1200,10 +1162,6 @@ class _InnerDerivativeJoiner(pymbolic.mapper.RecursiveMapper):
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):
@@ -1385,10 +1343,6 @@ class CollectorMixin(OperatorReducerMixin, LocalOpReducerMixin, FluxOpReducerMix
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
@@ -1433,12 +1387,9 @@ class BoundaryTagCollector(CollectorMixin, CombineMapper):
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
+# I'm not sure this works still.
+#class GeometricFactorCollector(CollectorMixin, CombineMapper):
+# pass
class BoundOperatorCollector(CSECachingMapperMixin, CollectorMixin, CombineMapper):
diff --git a/grudge/symbolic/operators.py b/grudge/symbolic/operators.py
index c74b083bbde5a2b9afefc638378099c624be2175..923092142e1d7ecf7ee11d596aced56d32c79054 100644
--- a/grudge/symbolic/operators.py
+++ b/grudge/symbolic/operators.py
@@ -35,6 +35,16 @@ from pytools import Record, memoize_method
# {{{ base classes
class Operator(pymbolic.primitives.Leaf):
+ """
+ .. attribute:: where
+
+ *None* for the default volume discretization or a boundary
+ tag for an operation on the denoted part of the boundary.
+ """
+
+ def __init__(self, where=None):
+ self.where = where
+
def stringifier(self):
from grudge.symbolic.mappers import StringifyMapper
return StringifyMapper
@@ -73,17 +83,15 @@ class Operator(pymbolic.primitives.Leaf):
return self.__class__ == other.__class__ and \
self.__getinitargs__() == other.__getinitargs__()
-
-class StatelessOperator(Operator):
def __getinitargs__(self):
- return ()
+ return (self.where,)
# }}}
# {{{ sum, integral, max
-class NodalReductionOperator(StatelessOperator):
+class NodalReductionOperator(Operator):
pass
@@ -159,12 +167,12 @@ class QuadratureStiffnessTOperator(DiffOperatorBase):
in favor of operators on the reference element.
"""
- def __init__(self, xyz_axis, quadrature_tag):
- DiffOperatorBase.__init__(self, xyz_axis)
- self.quadrature_tag = quadrature_tag
+ def __init__(self, xyz_axis, input_quadrature_tag, where=None):
+ super(QuadratureStiffnessTOperator, self).__init__(xyz_axis, where=where)
+ self.input_quadrature_tag = input_quadrature_tag
def __getinitargs__(self):
- return (self.xyz_axis, self.quadrature_tag)
+ return (self.xyz_axis, self.input_quadrature_tag, self.where)
mapper_method = intern("map_quad_stiffness_t")
@@ -179,16 +187,15 @@ def DiffOperatorVector(els):
# {{{ reference-element differentiation
class ReferenceDiffOperatorBase(Operator):
- def __init__(self, rst_axis):
- Operator.__init__(self)
+ def __init__(self, rst_axis, where=None):
+ super(ReferenceDiffOperatorBase, self).__init__(where)
self.rst_axis = rst_axis
+ if self.where is None:
+ 1/0
def __getinitargs__(self):
- return (self.rst_axis,)
-
- def preimage_ranges(self, eg):
- return eg.ranges
+ return (self.rst_axis, self.where)
def equal_except_for_axis(self, other):
return (type(self) == type(other)
@@ -197,18 +204,10 @@ class ReferenceDiffOperatorBase(Operator):
class ReferenceDifferentiationOperator(ReferenceDiffOperatorBase):
- @staticmethod
- def matrices(element_group):
- return element_group.differentiation_matrices
-
mapper_method = intern("map_ref_diff")
class ReferenceStiffnessTOperator(ReferenceDiffOperatorBase):
- @staticmethod
- def matrices(element_group):
- return element_group.stiffness_t_matrices
-
mapper_method = intern("map_ref_stiffness_t")
@@ -221,22 +220,17 @@ class ReferenceQuadratureStiffnessTOperator(ReferenceDiffOperatorBase):
when a :class:`StiffnessTOperator` is applied to a quadrature field.
"""
- def __init__(self, rst_axis, quadrature_tag):
- ReferenceDiffOperatorBase.__init__(self, rst_axis)
- self.quadrature_tag = quadrature_tag
+ def __init__(self, rst_axis, input_quadrature_tag, where=None):
+ ReferenceDiffOperatorBase.__init__(self, rst_axis, where)
+ self.input_quadrature_tag = input_quadrature_tag
def __getinitargs__(self):
- return (self.rst_axis, self.quadrature_tag)
+ return (self.rst_axis,
+ self.input_quadrature_tag,
+ self.where)
mapper_method = intern("map_ref_quad_stiffness_t")
- def preimage_ranges(self, eg):
- return eg.quadrature_info[self.quadrature_tag].ranges
-
- def matrices(self, element_group):
- return element_group.quadrature_info[self.quadrature_tag] \
- .ldis_quad_info.stiffness_t_matrices()
-
# }}}
# }}}
@@ -254,7 +248,7 @@ class ElementwiseLinearOperator(Operator):
mapper_method = intern("map_elementwise_linear")
-class ElementwiseMaxOperator(StatelessOperator):
+class ElementwiseMaxOperator(Operator):
mapper_method = intern("map_elementwise_max")
@@ -267,8 +261,9 @@ class QuadratureGridUpsampler(Operator):
In pre-processing, the boundary quad interpolation is specialized to
a separate operator, :class:`QuadratureBoundaryGridUpsampler`.
"""
- def __init__(self, quadrature_tag):
+ def __init__(self, quadrature_tag, where=None):
self.quadrature_tag = quadrature_tag
+ self.where = where
def __getinitargs__(self):
return (self.quadrature_tag,)
@@ -291,24 +286,6 @@ class QuadratureInteriorFacesGridUpsampler(Operator):
mapper_method = intern("map_quad_int_faces_grid_upsampler")
-
-class QuadratureBoundaryGridUpsampler(Operator):
- """
- .. note::
-
- This operator is purely for internal use. It is inserted
- by :class:`grudge.symbolic.mappers.OperatorSpecializer`
- when a :class:`MassOperator` is applied to a quadrature field.
- """
- def __init__(self, quadrature_tag, boundary_tag):
- self.quadrature_tag = quadrature_tag
- self.boundary_tag = boundary_tag
-
- def __getinitargs__(self):
- return (self.quadrature_tag, self.boundary_tag)
-
- mapper_method = intern("map_quad_bdry_grid_upsampler")
-
# }}}
@@ -345,7 +322,7 @@ class FilterOperator(ElementwiseLinearOperator):
return mat
-class OnesOperator(ElementwiseLinearOperator, StatelessOperator):
+class OnesOperator(ElementwiseLinearOperator):
def matrix(self, eg):
ldis = eg.local_discretization
@@ -353,7 +330,7 @@ class OnesOperator(ElementwiseLinearOperator, StatelessOperator):
return np.ones((node_count, node_count), dtype=np.float64)
-class AveragingOperator(ElementwiseLinearOperator, StatelessOperator):
+class AveragingOperator(ElementwiseLinearOperator):
def matrix(self, eg):
# average matrix, so that AVE*fields = cellaverage(fields)
# see Hesthaven and Warburton page 227
@@ -367,14 +344,14 @@ class AveragingOperator(ElementwiseLinearOperator, StatelessOperator):
return avg_mat
-class InverseVandermondeOperator(ElementwiseLinearOperator, StatelessOperator):
+class InverseVandermondeOperator(ElementwiseLinearOperator):
def matrix(self, eg):
return np.asarray(
la.inv(eg.local_discretization.vandermonde()),
order="C")
-class VandermondeOperator(ElementwiseLinearOperator, StatelessOperator):
+class VandermondeOperator(ElementwiseLinearOperator):
def matrix(self, eg):
return np.asarray(
eg.local_discretization.vandermonde(),
@@ -387,7 +364,7 @@ class VandermondeOperator(ElementwiseLinearOperator, StatelessOperator):
# {{{ mass operators
-class MassOperatorBase(ElementwiseLinearOperator, StatelessOperator):
+class MassOperatorBase(ElementwiseLinearOperator):
pass
@@ -445,11 +422,12 @@ class ReferenceQuadratureMassOperator(Operator):
when a :class:`MassOperator` is applied to a quadrature field.
"""
- def __init__(self, quadrature_tag):
+ def __init__(self, quadrature_tag, where=None):
+ super(Operator, self).__init__(self.where)
self.quadrature_tag = quadrature_tag
def __getinitargs__(self):
- return (self.quadrature_tag,)
+ return (self.quadrature_tag, self.where)
mapper_method = intern("map_ref_quad_mass")
diff --git a/grudge/symbolic/primitives.py b/grudge/symbolic/primitives.py
index 67543ed20a8eeae22a08fdfab5a8aa2a99775ed7..e5c8c2e276dde519c706e274a85f0d3edbb6980a 100644
--- a/grudge/symbolic/primitives.py
+++ b/grudge/symbolic/primitives.py
@@ -28,10 +28,13 @@ from six.moves import range
import numpy as np
import pymbolic.primitives
-from meshmode.mesh import BTAG_ALL
+from meshmode.mesh import BTAG_ALL, BTAG_NONE # noqa
from pymbolic.primitives import ( # noqa
- make_common_subexpression, If, Comparison)
+ cse_scope as cse_scope_base,
+ make_common_subexpression as cse, If, Comparison)
+from pymbolic.geometric_algebra import MultiVector
+from pytools.obj_array import join_fields, make_obj_array # noqa
class LeafBase(pymbolic.primitives.AlgebraicLeaf):
@@ -40,8 +43,48 @@ class LeafBase(pymbolic.primitives.AlgebraicLeaf):
return StringifyMapper
+class VTAG_ALL:
+ """This is used in a 'where' field to denote the volume discretization."""
+ pass
+
+
+__doc__ = """
+Symbols
+^^^^^^^
+
+.. autoclass:: Field
+.. autoclass:: make_sym_vector
+.. autoclass:: make_sym_array
+.. autoclass:: ScalarParameter
+.. autoclass:: CFunction
+
+Helpers
+^^^^^^^
+
+.. autoclass:: OperatorBinding
+.. autoclass:: PrioritizedSubexpression
+.. autoclass:: BoundaryPair
+.. autoclass:: Ones
+
+Geometry data
+^^^^^^^^^^^^^
+.. autoclass:: NodeCoordinateComponent
+.. autofunction:: nodes
+.. autofunction:: mv_nodes
+.. autofunction:: forward_metric_derivative
+.. autofunction:: inverse_metric_derivative
+.. autofunction:: area_element
+.. autofunction:: mv_normal
+.. autofunction:: normal
+"""
+
+
# {{{ variables
+class cse_scope(cse_scope_base): # noqa
+ DISCRETIZATION = "grudge_discretization"
+
+
Field = pymbolic.primitives.Variable
make_sym_vector = pymbolic.primitives.make_sym_vector
@@ -162,68 +205,64 @@ class BoundaryPair(LeafBase):
class Ones(LeafBase):
- def __init__(self, quadrature_tag=None):
+ def __init__(self, quadrature_tag=None, where=None):
+ self.where = where
self.quadrature_tag = quadrature_tag
def __getinitargs__(self):
- return (self.quadrature_tag,)
+ return (self.where, self.quadrature_tag,)
mapper_method = intern("map_ones")
# {{{ geometry data
-class NodeCoordinateComponent(LeafBase):
- def __init__(self, axis, quadrature_tag=None):
- self.axis = axis
- self.quadrature_tag = quadrature_tag
-
- def __getinitargs__(self):
- return (self.axis, self.quadrature_tag)
-
- mapper_method = intern("map_node_coordinate_component")
+class DiscretizationProperty(LeafBase):
+ """
+ .. attribute:: where
+ *None* for the default volume discretization or a boundary
+ tag for an operation on the denoted part of the boundary.
-def nodes(dim, quadrature_tag=None):
- return np.array([NodeCoordinateComponent(i, quadrature_tag)
- for i in range(dim)], dtype=object)
+ .. attribute:: quadrature_tag
+ quadrature tag for the grid on
+ which this geometric factor is needed, or None for
+ nodal representation.
+ """
-class BoundaryNormalComponent(LeafBase):
- def __init__(self, boundary_tag, axis, quadrature_tag=None):
- self.boundary_tag = boundary_tag
- self.axis = axis
+ def __init__(self, quadrature_tag, where=None):
self.quadrature_tag = quadrature_tag
+ self.where = where
def __getinitargs__(self):
- return (self.boundary_tag, self.axis, self.quadrature_tag)
+ return (self.quadrature_tag, self.where)
- mapper_method = intern("map_normal_component")
+class NodeCoordinateComponent(DiscretizationProperty):
-def normal(tag, dimensions):
- return np.array([BoundaryNormalComponent(tag, i)
- for i in range(dimensions)], dtype=object)
+ def __init__(self, axis, quadrature_tag=None, where=None):
+ super(NodeCoordinateComponent, self).__init__(quadrature_tag, where)
+ self.axis = axis
+
+ def __getinitargs__(self):
+ return (self.axis, self.quadrature_tag)
+ mapper_method = intern("map_node_coordinate_component")
-class GeometricFactorBase(LeafBase):
- def __init__(self, quadrature_tag):
- """
- :param quadrature_tag: quadrature tag for the grid on
- which this geometric factor is needed, or None for
- nodal representation.
- """
- self.quadrature_tag = quadrature_tag
- def __getinitargs__(self):
- return (self.quadrature_tag,)
+def nodes(ambient_dim, quadrature_tag=None, where=None):
+ return np.array([NodeCoordinateComponent(i, quadrature_tag, where)
+ for i in range(ambient_dim)], dtype=object)
-class Jacobian(GeometricFactorBase):
- mapper_method = intern("map_jacobian")
+def mv_nodes(ambient_dim, quadrature_tag=None, where=None):
+ return MultiVector(
+ nodes(ambient_dim, quadrature_tag=quadrature_tag, where=where))
-class ForwardMetricDerivative(GeometricFactorBase):
+def forward_metric_derivative(xyz_axis, rst_axis, where=None,
+ quadrature_tag=None):
r"""
Pointwise metric derivatives representing
@@ -231,48 +270,87 @@ class ForwardMetricDerivative(GeometricFactorBase):
\frac{d x_{\mathtt{xyz\_axis}} }{d r_{\mathtt{rst\_axis}} }
"""
+ from grudge.symbolic.operators import (
+ ReferenceDifferentiationOperator, QuadratureGridUpsampler)
+ diff_op = ReferenceDifferentiationOperator(
+ rst_axis, where=where)
- def __init__(self, quadrature_tag, xyz_axis, rst_axis):
- """
- :param quadrature_tag: quadrature tag for the grid on
- which this geometric factor is needed, or None for
- nodal representation.
- """
+ result = diff_op(NodeCoordinateComponent(xyz_axis, where=where))
- GeometricFactorBase.__init__(self, quadrature_tag)
- self.xyz_axis = xyz_axis
- self.rst_axis = rst_axis
+ if quadrature_tag is not None:
+ result = QuadratureGridUpsampler(quadrature_tag, where)(result)
- def __getinitargs__(self):
- return (self.quadrature_tag, self.xyz_axis, self.rst_axis)
+ return cse(
+ result,
+ prefix="dx%d_dr%d" % (xyz_axis, rst_axis),
+ scope=cse_scope.DISCRETIZATION)
- mapper_method = intern("map_forward_metric_derivative")
+def parametrization_derivative(ambient_dim, dim=None, where=None,
+ quadrature_tag=None):
+ if dim is None:
+ dim = ambient_dim
-class InverseMetricDerivative(GeometricFactorBase):
- r"""
- Pointwise metric derivatives representing
+ par_grad = np.zeros((ambient_dim, dim), np.object)
+ for i in range(ambient_dim):
+ for j in range(dim):
+ par_grad[i, j] = forward_metric_derivative(
+ i, j, where=where, quadrature_tag=quadrature_tag)
- .. math::
+ from pytools import product
+ return product(MultiVector(vec) for vec in par_grad.T)
- \frac{d r_{\mathtt{rst\_axis}} }{d x_{\mathtt{xyz\_axis}} }
- """
- def __init__(self, quadrature_tag, rst_axis, xyz_axis):
- """
- :param quadrature_tag: quadrature tag for the grid on
- which this geometric factor is needed, or None for
- nodal representation.
- """
+def inverse_metric_derivative(rst_axis, xyz_axis, ambient_dim, dim=None,
+ quadrature_tag=None):
+ if dim is None:
+ dim = ambient_dim
- GeometricFactorBase.__init__(self, quadrature_tag)
- self.rst_axis = rst_axis
- self.xyz_axis = xyz_axis
+ if dim != ambient_dim:
+ raise ValueError("not clear what inverse_metric_derivative means if "
+ "the derivative matrix is not square")
+ raise NotImplementedError()
+
+
+def pseudoscalar(ambient_dim, dim=None, where=None, quadrature_tag=None):
+ if dim is None:
+ dim = ambient_dim
+
+ return cse(
+ parametrization_derivative(ambient_dim, dim, where=where,
+ quadrature_tag=quadrature_tag)
+ .project_max_grade(),
+ "pseudoscalar", cse_scope.DISCRETIZATION)
+
+
+def area_element(ambient_dim, dim=None, where=None, quadrature_tag=None):
+ return cse(
+ CFunction("sqrt")(
+ pseudoscalar(ambient_dim, dim, where, quadrature_tag=quadrature_tag)
+ .norm_squared()),
+ "area_element", cse_scope.DISCRETIZATION)
+
+
+def mv_normal(tag, ambient_dim, dim=None, quadrature_tag=None):
+ """Exterior unit normal as a MultiVector."""
+
+ if dim is None:
+ dim = ambient_dim - 1
+
+ # Don't be tempted to add a sign here. As it is, it produces
+ # exterior normals for positively oriented curves.
+
+ pder = (
+ pseudoscalar(ambient_dim, dim, tag, quadrature_tag=quadrature_tag)
+ /
+ area_element(ambient_dim, dim, tag, quadrature_tag=quadrature_tag))
+ return cse(pder.I | pder, "normal",
+ cse_scope.DISCRETIZATION)
- def __getinitargs__(self):
- return (self.quadrature_tag, self.rst_axis, self.xyz_axis)
- mapper_method = intern("map_inverse_metric_derivative")
+def normal(tag, ambient_dim, dim=None, quadrature_tag=None):
+ return mv_normal(tag, ambient_dim, dim,
+ quadrature_tag=quadrature_tag).as_vector()
# }}}