From 79b0b082fed3587b76a6a11388d477dbea12974f Mon Sep 17 00:00:00 2001 From: Andreas Kloeckner Date: Mon, 20 Nov 2017 23:42:57 -0600 Subject: [PATCH] Rename grudge.Discretization to DGDiscretizationWithBoundaries, refactor sub-discr/conn getter interfaces --- examples/advection/var-velocity.py | 6 +- examples/advection/weak.py | 6 +- examples/wave/var-propagation-speed.py | 4 +- examples/wave/wave-min.py | 4 +- grudge/__init__.py | 2 +- grudge/discretization.py | 207 +++++++++++++++++++------ grudge/execution.py | 113 ++++---------- grudge/shortcuts.py | 15 +- grudge/symbolic/dofdesc_inference.py | 2 +- grudge/symbolic/primitives.py | 1 + test/test_grudge.py | 12 +- 11 files changed, 211 insertions(+), 161 deletions(-) diff --git a/examples/advection/var-velocity.py b/examples/advection/var-velocity.py index edce0760..ddc75352 100644 --- a/examples/advection/var-velocity.py +++ b/examples/advection/var-velocity.py @@ -27,7 +27,7 @@ from pyopencl.tools import ( # noqa import logging logger = logging.getLogger(__name__) -from grudge import sym, bind, Discretization +from grudge import sym, bind, DGDiscretizationWithBoundaries from pytools.obj_array import join_fields @@ -44,7 +44,7 @@ def main(write_output=True, order=4): dt_factor = 4 h = 1/20 - discr = Discretization(cl_ctx, mesh, order=order) + discr = DGDiscretizationWithBoundaries(cl_ctx, mesh, order=order) sym_x = sym.nodes(2) @@ -68,7 +68,7 @@ def main(write_output=True, order=4): from grudge.models.advection import VariableCoefficientAdvectionOperator - discr = Discretization(cl_ctx, mesh, order=order) + discr = DGDiscretizationWithBoundaries(cl_ctx, mesh, order=order) op = VariableCoefficientAdvectionOperator(2, advec_v, inflow_u=u_analytic(sym.nodes(dim, sym.BTAG_ALL)), flux_type=flux_type) diff --git a/examples/advection/weak.py b/examples/advection/weak.py index 967eb95f..cdadb49e 100644 --- a/examples/advection/weak.py +++ b/examples/advection/weak.py @@ -27,7 +27,7 @@ from pyopencl.tools import ( # noqa import logging logger = logging.getLogger(__name__) -from grudge import sym, bind, Discretization +from grudge import sym, bind, DGDiscretizationWithBoundaries import numpy.linalg as la @@ -48,7 +48,7 @@ def main(write_output=True, order=4): dt_factor = 4 h = 1/20 - discr = Discretization(cl_ctx, mesh, order=order) + discr = DGDiscretizationWithBoundaries(cl_ctx, mesh, order=order) c = np.array([0.1,0.1]) norm_c = la.norm(c) @@ -66,7 +66,7 @@ def main(write_output=True, order=4): from grudge.models.advection import WeakAdvectionOperator from meshmode.mesh import BTAG_ALL, BTAG_NONE - discr = Discretization(cl_ctx, mesh, order=order) + discr = DGDiscretizationWithBoundaries(cl_ctx, mesh, order=order) op = WeakAdvectionOperator(c, inflow_u=u_analytic(sym.nodes(dim, sym.BTAG_ALL)), flux_type=flux_type) diff --git a/examples/wave/var-propagation-speed.py b/examples/wave/var-propagation-speed.py index 6e1a7ea7..f161d034 100644 --- a/examples/wave/var-propagation-speed.py +++ b/examples/wave/var-propagation-speed.py @@ -28,7 +28,7 @@ THE SOFTWARE. import numpy as np import pyopencl as cl from grudge.shortcuts import set_up_rk4 -from grudge import sym, bind, Discretization +from grudge import sym, bind, DGDiscretizationWithBoundaries def main(write_output=True, order=4): @@ -42,7 +42,7 @@ def main(write_output=True, order=4): b=(0.5,)*dims, n=(20,)*dims) - discr = Discretization(cl_ctx, mesh, order=order) + discr = DGDiscretizationWithBoundaries(cl_ctx, mesh, order=order) source_center = np.array([0.1, 0.22, 0.33])[:mesh.dim] source_width = 0.05 diff --git a/examples/wave/wave-min.py b/examples/wave/wave-min.py index 6e2baa1b..89f4bfa8 100644 --- a/examples/wave/wave-min.py +++ b/examples/wave/wave-min.py @@ -28,7 +28,7 @@ THE SOFTWARE. import numpy as np import pyopencl as cl from grudge.shortcuts import set_up_rk4 -from grudge import sym, bind, Discretization +from grudge import sym, bind, DGDiscretizationWithBoundaries def main(write_output=True, order=4): @@ -49,7 +49,7 @@ def main(write_output=True, order=4): print("%d elements" % mesh.nelements) - discr = Discretization(cl_ctx, mesh, order=order) + discr = DGDiscretizationWithBoundaries(cl_ctx, mesh, order=order) source_center = np.array([0.1, 0.22, 0.33])[:mesh.dim] source_width = 0.05 diff --git a/grudge/__init__.py b/grudge/__init__.py index c6a5d6b7..b854007a 100644 --- a/grudge/__init__.py +++ b/grudge/__init__.py @@ -24,4 +24,4 @@ THE SOFTWARE. import grudge.sym # noqa from grudge.execution import bind # noqa -from grudge.discretization import Discretization # noqa +from grudge.discretization import DGDiscretizationWithBoundaries # noqa diff --git a/grudge/discretization.py b/grudge/discretization.py index b3596acb..cf640c2c 100644 --- a/grudge/discretization.py +++ b/grudge/discretization.py @@ -27,22 +27,15 @@ from pytools import memoize_method from grudge import sym +# FIXME Naming not ideal class DiscretizationBase(object): pass -class Discretization(DiscretizationBase): +class DGDiscretizationWithBoundaries(DiscretizationBase): """ - .. attribute :: volume_discr - - .. automethod :: boundary_connection - .. automethod :: boundary_discr - - .. automethod :: interior_faces_connection - .. automethod :: interior_faces_discr - - .. automethod :: all_faces_connection - .. automethod :: all_faces_discr + .. automethod :: discr_from_dd + .. automethod :: connection_from_dds .. autoattribute :: cl_context .. autoattribute :: dim @@ -67,8 +60,8 @@ class Discretization(DiscretizationBase): from meshmode.discretization import Discretization - self.volume_discr = Discretization(cl_ctx, mesh, - self.get_group_factory_for_quadrature_tag(sym.QTAG_NONE)) + self._volume_discr = Discretization(cl_ctx, mesh, + self.group_factory_for_quadrature_tag(sym.QTAG_NONE)) # {{{ management of discretization-scoped common subexpressions @@ -80,40 +73,147 @@ class Discretization(DiscretizationBase): # }}} - def get_group_factory_for_quadrature_tag(self, quadrature_tag): - if quadrature_tag is not sym.QTAG_NONE: - # FIXME - raise NotImplementedError("quadrature") + @memoize_method + def discr_from_dd(self, dd): + dd = sym.as_dofdesc(dd) + + qtag = dd.quadrature_tag + + if dd.is_volume(): + if qtag is not sym.QTAG_NONE: + # FIXME + raise NotImplementedError("quadrature") + return self._volume_discr + + elif dd.domain_tag is sym.FACE_RESTR_ALL: + return self._all_faces_discr(qtag) + elif dd.domain_tag is sym.FACE_RESTR_INTERIOR: + return self._interior_faces_discr(qtag) + elif dd.is_boundary(): + return self._boundary_discr(dd.domain_tag, qtag) + else: + raise ValueError("DOF desc tag not understood: " + str(dd)) + + @memoize_method + def connection_from_dds(self, from_dd, to_dd): + from_dd = sym.as_dofdesc(from_dd) + to_dd = sym.as_dofdesc(to_dd) + + if from_dd.quadrature_tag is not sym.QTAG_NONE: + raise ValueError("cannot interpolate *from* a quadrature grid") + + qtag = to_dd.quadrature_tag + + if from_dd.is_volume(): + if to_dd.domain_tag is sym.FACE_RESTR_ALL: + return self._all_faces_volume_connection(qtag) + if to_dd.domain_tag is sym.FACE_RESTR_INTERIOR: + return self._interior_faces_connection(qtag) + elif to_dd.is_boundary(): + if from_dd.quadrature_tag is sym.QTAG_NONE: + return self._boundary_connection(to_dd.domain_tag) + else: + from meshmode.connection import ChainedDiscretizationConnection + return ChainedDiscretizationConnection( + [self._boundary_connection(to_dd.domain_tag.tag), + self._boundary_to_quad_connection( + to_dd.domain_tag.tag, to_dd.quadrature_tag)]) + + else: + raise ValueError("cannot interpolate from volume to: " + str(to_dd)) + + elif from_dd.domain_tag is sym.FACE_RESTR_INTERIOR: + if to_dd.domain_tag is sym.FACE_RESTR_ALL and qtag is sym.QTAG_NONE: + return self._all_faces_connection(None) + else: + raise ValueError( + "cannot interpolate from interior faces to: " + + str(to_dd)) + + elif from_dd.is_boundary(): + if to_dd.domain_tag is sym.FACE_RESTR_ALL and qtag is sym.QTAG_NONE: + return self._all_faces_connection(from_dd.domain_tag) + else: + raise ValueError( + "cannot interpolate from interior faces to: " + + str(to_dd)) + + else: + raise ValueError("cannot interpolate from: " + str(from_dd)) + + def group_factory_for_quadrature_tag(self, quadrature_tag): + """ + OK to override in user code to control mode/node choice. + """ + + if quadrature_tag is None: + quadrature_tag = sym.QTAG_NONE from meshmode.discretization.poly_element import \ - PolynomialWarpAndBlendGroupFactory + PolynomialWarpAndBlendGroupFactory, \ + QuadratureSimplexGroupFactory - return PolynomialWarpAndBlendGroupFactory(order=self.order) + if quadrature_tag is not sym.QTAG_NONE: + return QuadratureSimplexGroupFactory(order=self.order) + else: + return PolynomialWarpAndBlendGroupFactory(order=self.order) + + @memoize_method + def _quad_volume_discr(self, quadrature_tag): + from meshmode.discretization import Discretization + + return Discretization(self._volume_discr.cl_context, self._volume_discr.mesh, + self.group_factory_for_quadrature_tag(quadrature_tag)) # {{{ boundary @memoize_method - def boundary_connection(self, boundary_tag, quadrature_tag=None): + def _boundary_connection(self, boundary_tag): from meshmode.discretization.connection import make_face_restriction return make_face_restriction( - self.volume_discr, - self.get_group_factory_for_quadrature_tag(quadrature_tag), + self._volume_discr, + self.group_factory_for_quadrature_tag(sym.QTAG_NONE), boundary_tag=boundary_tag) - def boundary_discr(self, boundary_tag, quadrature_tag=None): - return self.boundary_connection(boundary_tag, quadrature_tag).to_discr + @memoize_method + def _boundary_to_quad_connection(self, boundary_tag, quadrature_tag): + assert quadrature_tag is not None + assert quadrature_tag is not sym.QTAG_NONE + + from meshmode.discretization.connection.same_mesh import \ + make_same_mesh_connection + + return make_same_mesh_connection( + self._boundary_discr(boundary_tag, sym.QTAG_NONE), + self._boundary_discr(boundary_tag, quadrature_tag)) + + @memoize_method + def _boundary_discr(self, boundary_tag, quadrature_tag=None): + if quadrature_tag is None: + quadrature_tag = sym.QTAG_NONE + + if quadrature_tag is sym.QTAG_NONE: + return self._boundary_connection(boundary_tag).to_discr + else: + no_quad_bdry_discr = self.boundary_discr(boundary_tag, sym.QTAG_NONE) + + from meshmode.discretization import Discretization + return Discretization( + self._volume_discr.cl_context, + no_quad_bdry_discr.mesh, + self.group_factory_for_quadrature_tag(quadrature_tag)) # }}} # {{{ interior faces @memoize_method - def interior_faces_connection(self, quadrature_tag=None): + def _interior_faces_connection(self, quadrature_tag=None): from meshmode.discretization.connection import ( make_face_restriction, FACE_RESTR_INTERIOR) return make_face_restriction( - self.volume_discr, - self.get_group_factory_for_quadrature_tag(quadrature_tag), + self._volume_discr, + self.group_factory_for_quadrature_tag(quadrature_tag), FACE_RESTR_INTERIOR, # FIXME: This will need to change as soon as we support @@ -121,8 +221,8 @@ class Discretization(DiscretizationBase): # types. per_face_groups=False) - def interior_faces_discr(self, quadrature_tag=None): - return self.interior_faces_connection(quadrature_tag).to_discr + def _interior_faces_discr(self, quadrature_tag=None): + return self._interior_faces_connection(quadrature_tag).to_discr @memoize_method def opposite_face_connection(self, quadrature_tag): @@ -134,19 +234,19 @@ class Discretization(DiscretizationBase): make_opposite_face_connection return make_opposite_face_connection( - self.interior_faces_connection(quadrature_tag)) + self._interior_faces_connection(quadrature_tag)) # }}} # {{{ all-faces @memoize_method - def all_faces_volume_connection(self, quadrature_tag=None): + def _all_faces_volume_connection(self, quadrature_tag=None): from meshmode.discretization.connection import ( make_face_restriction, FACE_RESTR_ALL) return make_face_restriction( - self.volume_discr, - self.get_group_factory_for_quadrature_tag(quadrature_tag), + self._volume_discr, + self.group_factory_for_quadrature_tag(quadrature_tag), FACE_RESTR_ALL, # FIXME: This will need to change as soon as we support @@ -154,17 +254,17 @@ class Discretization(DiscretizationBase): # types. per_face_groups=False) - def all_faces_discr(self, quadrature_tag=None): - return self.all_faces_volume_connection(quadrature_tag).to_discr + def _all_faces_discr(self, quadrature_tag=None): + return self._all_faces_volume_connection(quadrature_tag).to_discr @memoize_method - def all_faces_connection(self, boundary_tag, quadrature_tag=None): + def _all_faces_connection(self, boundary_tag): """Return a :class:`meshmode.discretization.connection.DiscretizationConnection` that goes from either - :meth:`interior_faces_discr` (if *boundary_tag* is None) + :meth:`_interior_faces_discr` (if *boundary_tag* is None) or - :meth:`boundary_discr` (if *boundary_tag* is not None) + :meth:`_boundary_discr` (if *boundary_tag* is not None) to a discretization containing all the faces of the volume discretization. """ @@ -172,37 +272,36 @@ class Discretization(DiscretizationBase): make_face_to_all_faces_embedding if boundary_tag is None: - faces_conn = self.interior_faces_connection(quadrature_tag) + faces_conn = self._interior_faces_connection() else: - faces_conn = self.boundary_connection(boundary_tag, quadrature_tag) + faces_conn = self._boundary_connection(boundary_tag) - return make_face_to_all_faces_embedding( - faces_conn, self.all_faces_discr(quadrature_tag)) + return make_face_to_all_faces_embedding(faces_conn, self._all_faces_discr()) # }}} @property def cl_context(self): - return self.volume_discr.cl_context + return self._volume_discr.cl_context @property def dim(self): - return self.volume_discr.dim + return self._volume_discr.dim @property def ambient_dim(self): - return self.volume_discr.ambient_dim + return self._volume_discr.ambient_dim @property def mesh(self): - return self.volume_discr.mesh + return self._volume_discr.mesh def empty(self, queue=None, dtype=None, extra_dims=None, allocator=None): - return self.volume_discr.empty(queue, dtype, extra_dims=extra_dims, + return self._volume_discr.empty(queue, dtype, extra_dims=extra_dims, allocator=allocator) def zeros(self, queue, dtype=None, extra_dims=None, allocator=None): - return self.volume_discr.zeros(queue, dtype, extra_dims=extra_dims, + return self._volume_discr.zeros(queue, dtype, extra_dims=extra_dims, allocator=allocator) def is_volume_where(self, where): @@ -234,8 +333,16 @@ class PointsDiscretization(DiscretizationBase): def nodes(self): return self._nodes - @property - def volume_discr(self): + def discr_from_dd(self, dd): + dd = sym.as_dofdesc(dd) + + if dd.quadrature_tag is not sym.QTAG_NONE: + raise ValueError("quadrature discretization requested from " + "PointsDiscretization") + if dd.domain_tag is not sym.DTAG_VOLUME_ALL: + raise ValueError("non-volume discretization requested from " + "PointsDiscretization") + return self diff --git a/grudge/execution.py b/grudge/execution.py index 84fb24c5..99a3d4e5 100644 --- a/grudge/execution.py +++ b/grudge/execution.py @@ -43,45 +43,24 @@ class ExecutionMapper(mappers.Evaluator, mappers.LocalOpReducerMixin): def __init__(self, queue, context, bound_op): super(ExecutionMapper, self).__init__(context) - self.discr = bound_op.discr + self.discrwb = bound_op.discrwb self.bound_op = bound_op self.queue = queue - def get_discr(self, dd): - qtag = dd.quadrature_tag - if qtag is None: - # FIXME: Remove once proper quadrature support arrives - qtag = sym.QTAG_NONE - - if dd.is_volume(): - if qtag is not sym.QTAG_NONE: - # FIXME - raise NotImplementedError("quadrature") - return self.discr.volume_discr - - elif dd.domain_tag is sym.FACE_RESTR_ALL: - return self.discr.all_faces_discr(qtag) - elif dd.domain_tag is sym.FACE_RESTR_INTERIOR: - return self.discr.interior_faces_discr(qtag) - elif dd.is_boundary(): - return self.discr.boundary_discr(dd.domain_tag, qtag) - else: - raise ValueError("DOF desc tag not understood: " + str(dd)) - # {{{ expression mappings ------------------------------------------------- def map_ones(self, expr): if expr.dd.is_scalar(): return 1 - discr = self.get_discr(expr.dd) + discr = self.discrwb.discr_from_dd(expr.dd) result = discr.empty(self.queue, allocator=self.bound_op.allocator) result.fill(1) return result def map_node_coordinate_component(self, expr): - discr = self.get_discr(expr.dd) + discr = self.discrwb.discr_from_dd(expr.dd) return discr.nodes()[expr.axis].with_queue(self.queue) def map_grudge_variable(self, expr): @@ -89,7 +68,7 @@ class ExecutionMapper(mappers.Evaluator, value = self.context[expr.name] if not expr.dd.is_scalar() and isinstance(value, Number): - discr = self.get_discr(expr.dd) + discr = self.discrwb.discr_from_dd(expr.dd) ary = discr.empty(self.queue) ary.fill(value) value = ary @@ -104,7 +83,7 @@ class ExecutionMapper(mappers.Evaluator, from numbers import Number if not dd.is_scalar() and isinstance(value, Number): - discr = self.get_discr(dd) + discr = self.discrwb.discr_from_dd(dd) ary = discr.empty(self.queue) ary.fill(value) value = ary @@ -200,7 +179,7 @@ class ExecutionMapper(mappers.Evaluator, knl = lp.split_iname(knl, "i", 16, inner_tag="l.0") return lp.tag_inames(knl, dict(k="g.0")) - discr = self.get_discr(op.dd_in) + discr = self.discrwb.discr_from_dd(op.dd_in) # FIXME: This shouldn't really assume that it's dealing with a volume # input. What about quadrature? What about boundaries? @@ -237,46 +216,7 @@ class ExecutionMapper(mappers.Evaluator, return out def map_interpolation(self, op, field_expr): - if op.dd_in.quadrature_tag not in [None, sym.QTAG_NONE]: - raise ValueError("cannot interpolate *from* a quadrature grid") - - dd_in = op.dd_in - dd_out = op.dd_out - - qtag = dd_out.quadrature_tag - if qtag is None: - # FIXME: Remove once proper quadrature support arrives - qtag = sym.QTAG_NONE - - if dd_in.is_volume(): - if dd_out.domain_tag is sym.FACE_RESTR_ALL: - conn = self.discr.all_faces_connection(qtag) - elif dd_out.domain_tag is sym.FACE_RESTR_INTERIOR: - conn = self.discr.interior_faces_connection(qtag) - elif dd_out.is_boundary(): - conn = self.discr.boundary_connection(dd_out.domain_tag, qtag) - else: - raise ValueError("cannot interpolate from volume to: " + str(dd_out)) - - elif dd_in.domain_tag is sym.FACE_RESTR_INTERIOR: - if dd_out.domain_tag is sym.FACE_RESTR_ALL: - conn = self.discr.all_faces_connection(None, qtag) - else: - raise ValueError( - "cannot interpolate from interior faces to: " - + str(dd_out)) - - elif dd_in.is_boundary(): - if dd_out.domain_tag is sym.FACE_RESTR_ALL: - conn = self.discr.all_faces_connection(dd_in.domain_tag, qtag) - else: - raise ValueError( - "cannot interpolate from interior faces to: " - + str(dd_out)) - - else: - raise ValueError("cannot interpolate from: " + str(dd_in)) - + conn = self.discrwb.connection_from_dds(op.dd_in, op.dd_out) return conn(self.queue, self.rec(field_expr)).with_queue(self.queue) def map_opposite_interior_face_swap(self, op, field_expr): @@ -287,7 +227,7 @@ class ExecutionMapper(mappers.Evaluator, # FIXME: Remove once proper quadrature support arrives qtag = sym.QTAG_NONE - return self.discr.opposite_face_connection(qtag)( + return self.discrwb.opposite_face_connection(qtag)( self.queue, self.rec(field_expr)).with_queue(self.queue) # {{{ face mass operator @@ -313,12 +253,7 @@ class ExecutionMapper(mappers.Evaluator, knl = lp.split_iname(knl, "i", 16, inner_tag="l.0") return lp.tag_inames(knl, dict(k="g.0")) - qtag = op.dd_in.quadrature_tag - if qtag is None: - # FIXME: Remove once proper quadrature support arrives - qtag = sym.QTAG_NONE - - all_faces_conn = self.discr.all_faces_volume_connection(qtag) + all_faces_conn = self.discrwb.connection_from_dds("vol", op.dd_in) all_faces_discr = all_faces_conn.to_discr vol_discr = all_faces_conn.from_discr @@ -363,7 +298,7 @@ class ExecutionMapper(mappers.Evaluator, for name, expr in six.iteritems(kdescr.input_mappings): kwargs[name] = self.rec(expr) - discr = self.get_discr(kdescr.governing_dd) + discr = self.discrwb.discr_from_dd(kdescr.governing_dd) for name in kdescr.scalar_args(): v = kwargs[name] if isinstance(v, (int, float)): @@ -388,14 +323,14 @@ class ExecutionMapper(mappers.Evaluator, assignments = [] for name, expr in zip(insn.names, insn.exprs): value = self.rec(expr) - self.discr._discr_scoped_subexpr_name_to_value[name] = value + self.discrwb._discr_scoped_subexpr_name_to_value[name] = value assignments.append((name, value)) return assignments, [] def map_insn_assign_from_discr_scoped(self, insn): return [(insn.name, - self.discr._discr_scoped_subexpr_name_to_value[insn.name])], [] + self.discrwb._discr_scoped_subexpr_name_to_value[insn.name])], [] def map_insn_diff_batch_assign(self, insn): field = self.rec(insn.field) @@ -405,13 +340,13 @@ class ExecutionMapper(mappers.Evaluator, # This should be unified with map_elementwise_linear, which should # be extended to support batching. - discr = self.get_discr(repr_op.dd_in) + discr = self.discrwb.discr_from_dd(repr_op.dd_in) # FIXME: Enable # assert repr_op.dd_in == repr_op.dd_out assert repr_op.dd_in.domain_tag == repr_op.dd_out.domain_tag - @memoize_in(self.discr, "reference_derivative_knl") + @memoize_in(self.discrwb, "reference_derivative_knl") def knl(): knl = lp.make_kernel( """{[imatrix,k,i,j]: @@ -459,8 +394,8 @@ class ExecutionMapper(mappers.Evaluator, # {{{ bound operator class BoundOperator(object): - def __init__(self, discr, discr_code, eval_code, debug_flags, allocator=None): - self.discr = discr + def __init__(self, discrwb, discr_code, eval_code, debug_flags, allocator=None): + self.discrwb = discrwb self.discr_code = discr_code self.eval_code = eval_code self.operator_data_cache = {} @@ -490,13 +425,16 @@ class BoundOperator(object): from pytools.obj_array import with_object_array_or_scalar - # {{{ discr-scope evaluation + # {{{ discrwb-scope evaluation - if any(result_var.name not in self.discr._discr_scoped_subexpr_name_to_value + if any( + (result_var.name not in + self.discrwb._discr_scoped_subexpr_name_to_value) for result_var in self.discr_code.result): - # need to do discr-scope evaluation - discr_eval_context = {} - self.discr_code.execute(ExecutionMapper(queue, discr_eval_context, self)) + # need to do discrwb-scope evaluation + discrwb_eval_context = {} + self.discr_code.execute( + ExecutionMapper(queue, discrwb_eval_context, self)) # }}} @@ -504,7 +442,8 @@ class BoundOperator(object): for name, var in six.iteritems(context): new_context[name] = with_object_array_or_scalar(replace_queue, var) - return self.eval_code.execute(ExecutionMapper(queue, new_context, self)) + return self.eval_code.execute( + ExecutionMapper(queue, new_context, self)) # }}} diff --git a/grudge/shortcuts.py b/grudge/shortcuts.py index 0e6940ff..bb12689b 100644 --- a/grudge/shortcuts.py +++ b/grudge/shortcuts.py @@ -44,13 +44,16 @@ def set_up_rk4(field_var_name, dt, fields, rhs, t_start=0): return dt_stepper -def make_visualizer(discr, vis_order): +def make_visualizer(discrwb, 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) + with cl.CommandQueue(discrwb.cl_context) as queue: + return make_visualizer(queue, discrwb.discr_from_dd("vol"), vis_order) -def make_boundary_visualizer(discr, vis_order): +def make_boundary_visualizer(discrwb, 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) + from grudge import sym + with cl.CommandQueue(discrwb.cl_context) as queue: + return make_visualizer( + queue, discrwb.discr_from_dd(sym.BTAG_ALL), + vis_order) diff --git a/grudge/symbolic/dofdesc_inference.py b/grudge/symbolic/dofdesc_inference.py index 9cb54357..832c6a03 100644 --- a/grudge/symbolic/dofdesc_inference.py +++ b/grudge/symbolic/dofdesc_inference.py @@ -171,7 +171,7 @@ class DOFDescInferenceMapper(RecursiveMapper, CSECachingMapperMixin): " in '%s'" % ( type(expr).__name__, - op_dd, expr.dd_in, + op_dd, expr.op.dd_in, str(expr))) return operator.dd_out diff --git a/grudge/symbolic/primitives.py b/grudge/symbolic/primitives.py index 2f3e8787..d64aef3a 100644 --- a/grudge/symbolic/primitives.py +++ b/grudge/symbolic/primitives.py @@ -182,6 +182,7 @@ class DOFDesc(object): elif domain_tag is None: pass elif domain_tag in [BTAG_ALL, BTAG_REALLY_ALL, BTAG_NONE]: + # FIXME: Should wrap these in DTAG_BOUNDARY pass elif isinstance(domain_tag, DTAG_BOUNDARY): pass diff --git a/test/test_grudge.py b/test/test_grudge.py index bd47f823..66026990 100644 --- a/test/test_grudge.py +++ b/test/test_grudge.py @@ -37,7 +37,7 @@ from pyopencl.tools import ( # noqa import logging logger = logging.getLogger(__name__) -from grudge import sym, bind, Discretization +from grudge import sym, bind, DGDiscretizationWithBoundaries @pytest.mark.parametrize("dim", [2, 3]) @@ -64,7 +64,7 @@ def test_inverse_metric(ctx_factory, dim): from meshmode.mesh.processing import map_mesh mesh = map_mesh(mesh, m) - discr = Discretization(cl_ctx, mesh, order=4) + discr = DGDiscretizationWithBoundaries(cl_ctx, mesh, order=4) sym_op = ( sym.forward_metric_derivative_mat(mesh.dim) @@ -98,7 +98,7 @@ def test_1d_mass_mat_trig(ctx_factory): mesh = generate_regular_rect_mesh(a=(-4*np.pi,), b=(9*np.pi,), n=(17,), order=1) - discr = Discretization(cl_ctx, mesh, order=8) + discr = DGDiscretizationWithBoundaries(cl_ctx, mesh, order=8) x = sym.nodes(1) f = bind(discr, sym.cos(x[0])**2)(queue) @@ -139,7 +139,7 @@ def test_tri_diff_mat(ctx_factory, dim, order=4): mesh = generate_regular_rect_mesh(a=(-0.5,)*dim, b=(0.5,)*dim, n=(n,)*dim, order=4) - discr = Discretization(cl_ctx, mesh, order=4) + discr = DGDiscretizationWithBoundaries(cl_ctx, mesh, order=4) nabla = sym.nabla(dim) for axis in range(dim): @@ -182,7 +182,7 @@ def test_2d_gauss_theorem(ctx_factory): cl_ctx = cl.create_some_context() queue = cl.CommandQueue(cl_ctx) - discr = Discretization(cl_ctx, mesh, order=2) + discr = DGDiscretizationWithBoundaries(cl_ctx, mesh, order=2) def f(x): return sym.join_fields( @@ -272,7 +272,7 @@ def test_convergence_advec(ctx_factory, mesh_name, mesh_pars, op_type, flux_type from grudge.models.advection import ( StrongAdvectionOperator, WeakAdvectionOperator) - discr = Discretization(cl_ctx, mesh, order=order) + discr = DGDiscretizationWithBoundaries(cl_ctx, mesh, order=order) op_class = { "strong": StrongAdvectionOperator, "weak": WeakAdvectionOperator, -- GitLab