From b1eb9add18000cd1678603cda35689a3097170fa Mon Sep 17 00:00:00 2001
From: Andreas Kloeckner <inform@tiker.net>
Date: Tue, 5 Jan 2016 19:34:31 -0600
Subject: [PATCH] Redesign operator language to capture fluxes, include DOFDesc
---
examples/wave/wave-min.py | 2 +-
grudge/__init__.py | 3 -
grudge/discretization.py | 87 +++-
grudge/execution.py | 215 ++-------
grudge/models/wave.py | 83 ++--
grudge/symbolic/compiler.py | 282 +----------
grudge/symbolic/flux/__init__.py | 0
grudge/symbolic/flux/mappers.py | 188 --------
grudge/symbolic/flux/primitives.py | 256 ----------
grudge/symbolic/flux/tools.py | 85 ----
grudge/symbolic/mappers/__init__.py | 694 ++++++----------------------
grudge/symbolic/operators.py | 649 +++++++-------------------
grudge/symbolic/primitives.py | 476 +++++++++++++------
grudge/symbolic/tools.py | 56 +--
test/test_grudge.py | 8 +-
15 files changed, 834 insertions(+), 2250 deletions(-)
delete mode 100644 grudge/symbolic/flux/__init__.py
delete mode 100644 grudge/symbolic/flux/mappers.py
delete mode 100644 grudge/symbolic/flux/primitives.py
delete mode 100644 grudge/symbolic/flux/tools.py
diff --git a/examples/wave/wave-min.py b/examples/wave/wave-min.py
index 8fefcd6c..274b9b26 100644
--- a/examples/wave/wave-min.py
+++ b/examples/wave/wave-min.py
@@ -51,7 +51,7 @@ def main(write_output=True):
sym_source_center_dist = sym_x - source_center
sym_sin = sym.CFunction("sin")
sym_exp = sym.CFunction("sin")
- sym_t = sym.ScalarParameter("t")
+ sym_t = sym.ScalarVariable("t")
from grudge.models.wave import StrongWaveOperator
from meshmode.mesh import BTAG_ALL, BTAG_NONE
diff --git a/grudge/__init__.py b/grudge/__init__.py
index e156a2a9..c6a5d6b7 100644
--- a/grudge/__init__.py
+++ b/grudge/__init__.py
@@ -23,8 +23,5 @@ THE SOFTWARE.
"""
import grudge.sym # noqa
-import grudge.symbolic.flux.primitives as sym_flux # noqa
-
from grudge.execution import bind # noqa
-
from grudge.discretization import Discretization # noqa
diff --git a/grudge/discretization.py b/grudge/discretization.py
index fe847a3c..c992e471 100644
--- a/grudge/discretization.py
+++ b/grudge/discretization.py
@@ -24,6 +24,7 @@ THE SOFTWARE.
from pytools import memoize_method
+from grudge import sym
class Discretization(object):
@@ -36,6 +37,9 @@ class Discretization(object):
.. automethod :: interior_faces_connection
.. automethod :: interior_faces_discr
+ .. automethod :: all_faces_connection
+ .. automethod :: all_faces_discr
+
.. autoattribute :: cl_context
.. autoattribute :: dim
.. autoattribute :: ambient_dim
@@ -64,44 +68,101 @@ class Discretization(object):
self.order = order
self.quad_min_degrees = quad_min_degrees
+ # {{{ boundary
+
@memoize_method
def boundary_connection(self, boundary_tag, quadrature_tag=None):
from meshmode.discretization.poly_element import \
PolynomialWarpAndBlendGroupFactory
- if quadrature_tag is not None:
+ if quadrature_tag is not sym.QTAG_NONE:
+ # FIXME
raise NotImplementedError("quadrature")
from meshmode.discretization.connection import make_face_restriction
- _, _, boundary_connection = \
- make_face_restriction(
+ return make_face_restriction(
self.volume_discr,
PolynomialWarpAndBlendGroupFactory(order=self.order),
boundary_tag=boundary_tag)
- return self.bounds
-
def boundary_discr(self, boundary_tag, quadrature_tag=None):
return self.boundary_connection(boundary_tag, quadrature_tag).to_discr
+ # }}}
+
+ # {{{ interior faces
+
@memoize_method
def interior_faces_connection(self, quadrature_tag=None):
from meshmode.discretization.poly_element import \
PolynomialWarpAndBlendGroupFactory
- if quadrature_tag is not None:
+ if quadrature_tag is not sym.QTAG_NONE:
+ # FIXME
raise NotImplementedError("quadrature")
- from meshmode.discretization.connection import make_face_restriction
- _, _, boundary_connection = \
- make_face_restriction(
+ from meshmode.discretization.connection import (
+ make_face_restriction, FRESTR_INTERIOR_FACES)
+ return make_face_restriction(
self.volume_discr,
- PolynomialWarpAndBlendGroupFactory(order=self.order))
+ PolynomialWarpAndBlendGroupFactory(order=self.order),
+ FRESTR_INTERIOR_FACES,
- return self.bounds
+ # FIXME: This will need to change as soon as we support
+ # pyramids or other elements with non-identical face
+ # types.
+ per_face_groups=False)
- def interior_faces_discr(self, boundary_tag, quadrature_tag=None):
- return self.boundary_connection(boundary_tag, quadrature_tag).to_discr
+ def interior_faces_discr(self, quadrature_tag=None):
+ return self.interior_faces_connection(quadrature_tag).to_discr
+
+ # }}}
+
+ # {{{ all-faces
+
+ @memoize_method
+ def all_faces_discr(self, quadrature_tag=None):
+ if quadrature_tag is not sym.QTAG_NONE:
+ # FIXME
+ raise NotImplementedError("quadrature")
+
+ from meshmode.discretization.poly_element import \
+ PolynomialWarpAndBlendGroupFactory
+ from meshmode.discretization.connection import (
+ make_face_restriction, FRESTR_ALL_FACES)
+ return make_face_restriction(
+ self.volume_discr,
+ PolynomialWarpAndBlendGroupFactory(order=self.order),
+ FRESTR_ALL_FACES,
+
+ # FIXME: This will need to change as soon as we support
+ # pyramids or other elements with non-identical face
+ # types.
+ per_face_groups=False)
+
+ @memoize_method
+ def all_faces_connection(self, boundary_tag, quadrature_tag=None):
+ """Return a
+ :class:`meshmode.discretization.connection.DiscretizationConnection`
+ that goes from either
+ :meth:`interior_faces_discr` (if *boundary_tag* is None)
+ or
+ :meth:`boundary_discr` (if *boundary_tag* is not None)
+ to a discretization containing all the faces of the volume
+ discretization.
+ """
+ from meshmode.discretization.connection import \
+ make_face_to_all_faces_embedding
+
+ if boundary_tag is None:
+ faces_conn = self.interior_faces_connection(quadrature_tag)
+ else:
+ faces_conn = self.boundary_connection(boundary_tag, quadrature_tag)
+
+ return make_face_to_all_faces_embedding(
+ faces_conn, self.all_faces_discr(quadrature_tag))
+
+ # }}}
@property
def cl_context(self):
diff --git a/grudge/execution.py b/grudge/execution.py
index 2cf24b31..7531698f 100644
--- a/grudge/execution.py
+++ b/grudge/execution.py
@@ -47,27 +47,38 @@ class ExecutionMapper(mappers.Evaluator,
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.FRESTR_ALL_FACES:
+ return self.discr.all_faces_discr(qtag)
+ elif dd.domain_tag is sym.FRESTR_INTERIOR_FACES:
+ return self.discr.interior_faces_discr(qtag)
+ elif dd.is_boundary():
+ return self.discr.boundary_connection(dd.domain_tag, qtag)
+ else:
+ raise ValueError("DOF desc tag not understood: " + str(dd))
+
# {{{ expression mappings -------------------------------------------------
def map_ones(self, expr):
- if expr.quadrature_tag is not None:
- # FIXME
- raise NotImplementedError("ones on quad. grids")
+ discr = self.get_discr(expr.dd)
- result = self.discr.empty(self.queue, allocator=self.bound_op.allocator)
+ result = discr.empty(self.queue, allocator=self.bound_op.allocator)
result.fill(1)
return result
def map_node_coordinate_component(self, expr):
- if expr.quadrature_tag is not None:
- # FIXME
- raise NotImplementedError("node coordinate components on quad. grids")
-
- if self.discr.is_volume_where(expr.where):
- discr = self.discr.volume_discr
- else:
- discr = self.discr.boundary_discr
-
+ discr = self.get_discr(expr.dd)
return discr.nodes()[expr.axis].with_queue(self.queue)
def map_boundarize(self, op, field_expr):
@@ -75,7 +86,7 @@ class ExecutionMapper(mappers.Evaluator,
self.rec(field_expr), tag=op.tag,
kind=self.discr.compute_kind)
- def map_scalar_parameter(self, expr):
+ def map_grudge_variable(self, expr):
return self.context[expr.name]
def map_call(self, expr):
@@ -127,7 +138,7 @@ class ExecutionMapper(mappers.Evaluator,
"""{[k,i,j]:
0<=k<nelements and
0<=i,j<ndiscr_nodes}""",
- "result[k,i] = sum(j, diff_mat[i, j] * vec[k, j])",
+ "result[k,i] = sum(j, mat[i, j] * vec[k, j])",
default_offset=lp.auto, name="diff")
knl = lp.split_iname(knl, "i", 16, inner_tag="l.0")
@@ -154,7 +165,7 @@ class ExecutionMapper(mappers.Evaluator,
self.bound_op.elwise_linear_cache[cache_key] = matrix
- knl()(self.queue, diff_mat=matrix, result=grp.view(result),
+ knl()(self.queue, mat=matrix, result=grp.view(result),
vec=grp.view(field))
return result
@@ -283,134 +294,10 @@ class ExecutionMapper(mappers.Evaluator,
return zip(compiled.result_names(),
compiled(self, stats_callback)), []
- def exec_flux_batch_assign(self, insn):
- from pymbolic.primitives import is_zero
-
- class ZeroSpec:
- pass
-
- class BoundaryZeros(ZeroSpec):
- pass
-
- class VolumeZeros(ZeroSpec):
- pass
-
- def eval_arg(arg_spec):
- arg_expr, is_int = arg_spec
- arg = self.rec(arg_expr)
- if is_zero(arg):
- if insn.is_boundary and not is_int:
- return BoundaryZeros()
- else:
- return VolumeZeros()
- else:
- return arg
-
- args = [eval_arg(arg_expr)
- for arg_expr in insn.flux_var_info.arg_specs]
-
- from pytools import common_dtype
- max_dtype = common_dtype(
- [a.dtype for a in args if not isinstance(a, ZeroSpec)],
- self.discr.default_scalar_type)
-
- def cast_arg(arg):
- if isinstance(arg, BoundaryZeros):
- return self.discr.boundary_zeros(
- insn.repr_op.boundary_tag, dtype=max_dtype)
- elif isinstance(arg, VolumeZeros):
- return self.discr.volume_zeros(
- dtype=max_dtype)
- elif isinstance(arg, np.ndarray):
- return np.asarray(arg, dtype=max_dtype)
- else:
- return arg
-
- args = [cast_arg(arg) for arg in args]
-
- if insn.quadrature_tag is None:
- if insn.is_boundary:
- face_groups = self.discr.get_boundary(insn.repr_op.boundary_tag)\
- .face_groups
- else:
- face_groups = self.discr.face_groups
- else:
- if insn.is_boundary:
- face_groups = self.discr.get_boundary(insn.repr_op.boundary_tag)\
- .get_quadrature_info(insn.quadrature_tag).face_groups
- else:
- face_groups = self.discr.get_quadrature_info(insn.quadrature_tag) \
- .face_groups
-
- result = []
-
- for fg in face_groups:
- # grab module
- module = insn.get_module(self.discr, max_dtype)
- func = module.gather_flux
-
- # set up argument structure
- arg_struct = module.ArgStruct()
- for arg_name, arg in zip(insn.flux_var_info.arg_names, args):
- setattr(arg_struct, arg_name, arg)
- for arg_num, scalar_arg_expr in enumerate(
- insn.flux_var_info.scalar_parameters):
- setattr(arg_struct,
- "_scalar_arg_%d" % arg_num,
- self.rec(scalar_arg_expr))
-
- fof_shape = (fg.face_count*fg.face_length()*fg.element_count(),)
- all_fluxes_on_faces = [
- np.zeros(fof_shape, dtype=max_dtype)
- for f in insn.expressions]
- for i, fof in enumerate(all_fluxes_on_faces):
- setattr(arg_struct, "flux%d_on_faces" % i, fof)
-
- # make sure everything ended up in Boost.Python attributes
- # (i.e. empty __dict__)
- assert not arg_struct.__dict__, arg_struct.__dict__.keys()
-
- # perform gather
- func(fg, arg_struct)
-
- # do lift, produce output
- for name, flux_bdg, fluxes_on_faces in zip(insn.names, insn.expressions,
- all_fluxes_on_faces):
-
- if insn.quadrature_tag is None:
- if flux_bdg.op.is_lift:
- mat = fg.ldis_loc.lifting_matrix()
- scaling = fg.local_el_inverse_jacobians
- else:
- mat = fg.ldis_loc.multi_face_mass_matrix()
- scaling = None
- else:
- assert not flux_bdg.op.is_lift
- mat = fg.ldis_loc_quad_info.multi_face_mass_matrix()
- scaling = None
-
- out = self.discr.volume_zeros(dtype=fluxes_on_faces.dtype)
- self.bound_op.lift_flux(fg, mat, scaling, fluxes_on_faces, out)
-
- if self.discr.instrumented:
- from grudge.tools import lift_flops
-
- # correct for quadrature, too.
- self.discr.lift_flop_counter.add(lift_flops(fg))
-
- result.append((name, out))
-
- if not face_groups:
- # No face groups? Still assign context variables.
- for name, flux_bdg in zip(insn.names, insn.expressions):
- result.append((name, self.discr.volume_zeros()))
-
- return result, []
-
def exec_diff_batch_assign(self, insn):
field = self.rec(insn.field)
repr_op = insn.operators[0]
- if not isinstance(repr_op, sym.ReferenceDifferentiationOperator):
+ if not isinstance(repr_op, sym.RefDiffOperator):
# FIXME
raise NotImplementedError()
@@ -418,10 +305,8 @@ class ExecutionMapper(mappers.Evaluator,
# execution-wise.
# This should be unified with map_elementwise_linear, which should
# be extended to support batching.
- if self.discr.is_volume_where(repr_op.where):
- discr = self.discr.volume_discr
- else:
- discr = self.discr.boundary_discr
+
+ discr = self.get_discr(repr_op.dd_in)
return [
(name, discr.num_reference_derivative(
@@ -429,8 +314,6 @@ class ExecutionMapper(mappers.Evaluator,
.with_queue(self.queue))
for name, op in zip(insn.names, insn.operators)], []
- exec_quad_diff_batch_assign = exec_diff_batch_assign
-
# }}}
# }}}
@@ -469,12 +352,9 @@ class BoundOperator(object):
# {{{ process_sym_operator function
def process_sym_operator(sym_operator, post_bind_mapper=None,
- dumper=lambda name, sym_operator: None, mesh=None,
- type_hints={}):
+ dumper=lambda name, sym_operator: None, mesh=None):
import grudge.symbolic.mappers as mappers
- from grudge.symbolic.mappers.bc_to_flux import BCToFluxRewriter
- from grudge.symbolic.mappers.type_inference import TypeInferrer
dumper("before-bind", sym_operator)
sym_operator = mappers.OperatorBinder()(sym_operator)
@@ -492,22 +372,19 @@ def process_sym_operator(sym_operator, post_bind_mapper=None,
dumper("before-cfold", sym_operator)
sym_operator = mappers.CommutativeConstantFoldingMapper()(sym_operator)
- dumper("before-bc2flux", sym_operator)
- sym_operator = BCToFluxRewriter()(sym_operator)
-
# Ordering restriction:
#
# - Must run constant fold before first type inference pass, because zeros,
# while allowed, violate typing constraints (because they can't be assigned
# a unique type), and need to be killed before the type inferrer sees them.
- #
- # - Must run BC-to-flux before first type inferrer run so that zeros in
- # flux arguments can be removed.
- dumper("before-specializer", sym_operator)
- sym_operator = mappers.OperatorSpecializer(
- TypeInferrer()(sym_operator, type_hints)
- )(sym_operator)
+ # FIXME: Reenable type inference
+
+ # from grudge.symbolic.mappers.type_inference import TypeInferrer
+ # dumper("before-specializer", sym_operator)
+ # sym_operator = mappers.OperatorSpecializer(
+ # TypeInferrer()(sym_operator)
+ # )(sym_operator)
# Ordering restriction:
#
@@ -531,11 +408,9 @@ def process_sym_operator(sym_operator, post_bind_mapper=None,
dumper("before-cfold-2", sym_operator)
sym_operator = mappers.CommutativeConstantFoldingMapper()(sym_operator)
- dumper("before-derivative-join", sym_operator)
- sym_operator = mappers.DerivativeJoiner()(sym_operator)
-
- dumper("before-boundary-combiner", sym_operator)
- sym_operator = mappers.BoundaryCombiner(mesh)(sym_operator)
+ # FIXME: Reenable derivative joiner
+ # dumper("before-derivative-join", sym_operator)
+ # sym_operator = mappers.DerivativeJoiner()(sym_operator)
dumper("process-finished", sym_operator)
@@ -544,7 +419,7 @@ def process_sym_operator(sym_operator, post_bind_mapper=None,
# }}}
-def bind(discr, sym_operator, post_bind_mapper=lambda x: x, type_hints={},
+def bind(discr, sym_operator, post_bind_mapper=lambda x: x,
debug_flags=set(), allocator=None):
# from grudge.symbolic.mappers import QuadratureUpsamplerRemover
# sym_operator = QuadratureUpsamplerRemover(self.quad_min_degrees)(
@@ -563,12 +438,10 @@ def bind(discr, sym_operator, post_bind_mapper=lambda x: x, type_hints={},
sym_operator = process_sym_operator(sym_operator,
post_bind_mapper=post_bind_mapper,
dumper=dump_optemplate,
- mesh=discr.mesh,
- type_hints=type_hints)
+ mesh=discr.mesh)
from grudge.symbolic.compiler import OperatorCompiler
- code = OperatorCompiler()(sym_operator,
- type_hints=type_hints)
+ code = OperatorCompiler()(sym_operator)
bound_op = BoundOperator(discr, code,
debug_flags=debug_flags, allocator=allocator)
diff --git a/grudge/models/wave.py b/grudge/models/wave.py
index 9591573e..8d559426 100644
--- a/grudge/models/wave.py
+++ b/grudge/models/wave.py
@@ -31,7 +31,7 @@ import numpy as np
from grudge.models import HyperbolicOperator
from grudge.models.second_order import CentralSecondDerivative
from meshmode.mesh import BTAG_ALL, BTAG_NONE
-from grudge import sym, sym_flux
+from grudge import sym
from pytools.obj_array import join_fields
@@ -55,16 +55,16 @@ class StrongWaveOperator(HyperbolicOperator):
:math:`c` is assumed to be constant across all space.
"""
- def __init__(self, c, dimensions, source_f=0,
+ def __init__(self, c, ambient_dim, source_f=0,
flux_type="upwind",
dirichlet_tag=BTAG_ALL,
dirichlet_bc_f=0,
neumann_tag=BTAG_NONE,
radiation_tag=BTAG_NONE):
- assert isinstance(dimensions, int)
+ assert isinstance(ambient_dim, int)
self.c = c
- self.dimensions = dimensions
+ self.ambient_dim = ambient_dim
self.source_f = source_f
if self.c > 0:
@@ -80,12 +80,10 @@ class StrongWaveOperator(HyperbolicOperator):
self.flux_type = flux_type
- def flux(self):
- dim = self.dimensions
- w = sym_flux.FluxVectorPlaceholder(1+dim)
+ def flux(self, w):
u = w[0]
v = w[1:]
- normal = sym_flux.normal(dim)
+ normal = sym.normal(w.dd, self.ambient_dim)
flux_weak = join_fields(
np.dot(v.avg, normal),
@@ -108,17 +106,17 @@ class StrongWaveOperator(HyperbolicOperator):
return -self.c*flux_strong
def sym_operator(self):
- d = self.dimensions
+ d = self.ambient_dim
- w = sym.make_sym_vector("w", d+1)
+ w = sym.make_sym_array("w", d+1)
u = w[0]
v = w[1:]
# boundary conditions -------------------------------------------------
# dirichlet BCs -------------------------------------------------------
- dir_u = sym.RestrictToBoundary(self.dirichlet_tag) * u
- dir_v = sym.RestrictToBoundary(self.dirichlet_tag) * v
+ dir_u = sym.cse(sym.interp("vol", self.dirichlet_tag)(u))
+ dir_v = sym.cse(sym.interp("vol", self.dirichlet_tag)(v))
if self.dirichlet_bc_f:
# FIXME
from warnings import warn
@@ -130,25 +128,30 @@ class StrongWaveOperator(HyperbolicOperator):
else:
dir_bc = join_fields(-dir_u, dir_v)
+ dir_bc = sym.cse(dir_bc, "dir_bc")
+
# neumann BCs ---------------------------------------------------------
- neu_u = sym.RestrictToBoundary(self.neumann_tag) * u
- neu_v = sym.RestrictToBoundary(self.neumann_tag) * v
- neu_bc = join_fields(neu_u, -neu_v)
+ neu_u = sym.cse(sym.interp("vol", self.neumann_tag)(u))
+ neu_v = sym.cse(sym.interp("vol", self.neumann_tag)(v))
+ neu_bc = sym.cse(join_fields(neu_u, -neu_v), "neu_bc")
# radiation BCs -------------------------------------------------------
rad_normal = sym.normal(self.radiation_tag, d)
- rad_u = sym.RestrictToBoundary(self.radiation_tag) * u
- rad_v = sym.RestrictToBoundary(self.radiation_tag) * v
+ rad_u = sym.cse(sym.interp("vol", self.radiation_tag)(u))
+ rad_v = sym.cse(sym.interp("vol", self.radiation_tag)(v))
- rad_bc = join_fields(
+ rad_bc = sym.cse(join_fields(
0.5*(rad_u - self.sign*np.dot(rad_normal, rad_v)),
0.5*rad_normal*(np.dot(rad_normal, rad_v) - self.sign*rad_u)
- )
+ ), "rad_bc")
# entire operator -----------------------------------------------------
nabla = sym.nabla(d)
- flux_op = sym.get_flux_operator(self.flux())
+
+ def flux(pair):
+ return sym.interp(pair.dd, "all_faces")(
+ self.flux(pair))
result = (
- join_fields(
@@ -156,12 +159,13 @@ class StrongWaveOperator(HyperbolicOperator):
-self.c*(nabla*u)
)
+
- sym.InverseMassOperator() * (
- flux_op(w)
- + flux_op(sym.BoundaryPair(w, dir_bc, self.dirichlet_tag))
- + flux_op(sym.BoundaryPair(w, neu_bc, self.neumann_tag))
- + flux_op(sym.BoundaryPair(w, rad_bc, self.radiation_tag))
- ))
+ sym.InverseMassOperator()(
+ sym.FaceMassOperator()(
+ flux(sym.int_fpair(w))
+ + flux(sym.bv_fpair(self.dirichlet_tag, w, dir_bc))
+ + flux(sym.bv_fpair(self.neumann_tag, w, neu_bc))
+ + flux(sym.bv_fpair(self.radiation_tag, w, rad_bc))
+ )))
result[0] += self.source_f
@@ -196,7 +200,7 @@ class VariableVelocityStrongWaveOperator(HyperbolicOperator):
"""
def __init__(
- self, c, dimensions, source=0,
+ self, c, ambient_dim, source=0,
flux_type="upwind",
dirichlet_tag=BTAG_ALL,
neumann_tag=BTAG_NONE,
@@ -207,11 +211,11 @@ class VariableVelocityStrongWaveOperator(HyperbolicOperator):
):
"""*c* and *source* are optemplate expressions.
"""
- assert isinstance(dimensions, int)
+ assert isinstance(ambient_dim, int)
self.c = c
self.time_sign = time_sign
- self.dimensions = dimensions
+ self.ambient_dim = ambient_dim
self.source = source
self.dirichlet_tag = dirichlet_tag
@@ -224,13 +228,12 @@ class VariableVelocityStrongWaveOperator(HyperbolicOperator):
self.diffusion_scheme = diffusion_scheme
# {{{ flux ----------------------------------------------------------------
- def flux(self):
- dim = self.dimensions
- w = sym_flux.FluxVectorPlaceholder(2+dim)
- c = w[0]
- u = w[1]
- v = w[2:]
- normal = sym_flux.normal(dim)
+ def flux(self, w, c_vol):
+ u = w[0]
+ v = w[1:]
+ normal = sym.normal(w.tag, self.ambient_dim)
+
+ c = sym.RestrictToBoundary(w.tag)(c_vol)
flux = self.time_sign*1/2*join_fields(
c.ext * np.dot(v.ext, normal)
@@ -263,9 +266,9 @@ class VariableVelocityStrongWaveOperator(HyperbolicOperator):
return do
def sym_operator(self, with_sensor=False):
- d = self.dimensions
+ d = self.ambient_dim
- w = sym.make_sym_vector("w", d+1)
+ w = sym.make_sym_array("w", d+1)
u = w[0]
v = w[1:]
@@ -320,14 +323,14 @@ class VariableVelocityStrongWaveOperator(HyperbolicOperator):
# strong_form here allows the reuse the value of grad u.
grad_tgt = SecondDerivativeTarget(
- self.dimensions, strong_form=True,
+ self.ambient_dim, strong_form=True,
operand=arg)
self.diffusion_scheme.grad(grad_tgt, bc_getter=None,
dirichlet_tags=[], neumann_tags=[])
div_tgt = SecondDerivativeTarget(
- self.dimensions, strong_form=False,
+ self.ambient_dim, strong_form=False,
operand=diffusion_coeff*grad_tgt.minv_all)
self.diffusion_scheme.div(div_tgt,
diff --git a/grudge/symbolic/compiler.py b/grudge/symbolic/compiler.py
index 55c5fa6d..37dba746 100644
--- a/grudge/symbolic/compiler.py
+++ b/grudge/symbolic/compiler.py
@@ -25,7 +25,7 @@ THE SOFTWARE.
"""
-import six
+import six # noqa
from six.moves import zip, reduce
from pytools import Record, memoize_method, memoize
from grudge import sym
@@ -140,67 +140,6 @@ class Assign(Instruction):
return exec_mapper.exec_assign
-class FluxBatchAssign(Instruction):
- """
- .. attribute:: names
- .. attribute:: expressions
-
- A list of :class:`grudge.symbolic.primitives.OperatorBinding`
- instances bound to flux operators.
-
- .. note ::
-
- All operators in :attr:`expressions` are guaranteed to
- yield the same operator from
- :meth:`grudge.symbolic.operators.FluxOperatorBase.repr_op`.
-
- .. attribute:: repr_op
-
- The *repr_op* on which all operators agree.
-
- .. attribute:: quadrature_tag
- """
-
- def get_assignees(self):
- return set(self.names)
-
- @memoize_method
- def get_dependencies(self):
- deps = set()
- for wdflux in self.expressions:
- deps |= set(wdflux.interior_deps)
- deps |= set(wdflux.boundary_deps)
-
- dep_mapper = _make_dep_mapper(include_subscripts=False)
-
- from pytools import flatten
- return set(flatten(dep_mapper(dep) for dep in deps))
-
- def __str__(self):
- from grudge.symbolic.flux.mappers import PrettyFluxStringifyMapper as PFSM
- flux_strifier = PFSM()
- op_strifier = mappers.StringifyMapper(flux_stringify_mapper=flux_strifier)
-
- from pymbolic.mapper.stringifier import PREC_NONE
-
- lines = []
- lines.append("{ /* %s */" % self.repr_op)
-
- lines_expr = []
- for n, f in zip(self.names, self.expressions):
- lines_expr.append(" %s <- %s" % (n, op_strifier(f, PREC_NONE)))
-
- for n, str_f in getattr(flux_strifier, "cse_name_list", []):
- lines.append(" (flux-local) %s <- %s" % (n, str_f))
-
- lines.extend(lines_expr)
- lines.append("}")
- return "\n".join(lines)
-
- def get_execution_method(self, exec_mapper):
- return exec_mapper.exec_flux_batch_assign
-
-
class DiffBatchAssign(Instruction):
"""
:ivar names:
@@ -240,11 +179,6 @@ class DiffBatchAssign(Instruction):
return exec_mapper.exec_diff_batch_assign
-class QuadratureDiffBatchAssign(DiffBatchAssign):
- def get_execution_method(self, exec_mapper):
- return exec_mapper.exec_quad_diff_batch_assign
-
-
class FluxExchangeBatchAssign(Instruction):
"""
.. attribute:: names
@@ -605,12 +539,6 @@ class Code(object):
# {{{ compiler
class OperatorCompiler(mappers.IdentityMapper):
- class FluxRecord(Record):
- __slots__ = ["flux_expr", "dependencies", "repr_op"]
-
- class FluxBatch(Record):
- __slots__ = ["flux_exprs", "repr_op"]
-
def __init__(self, prefix="_expr", max_vectors_in_batch_expr=None):
super(OperatorCompiler, self).__init__()
self.prefix = prefix
@@ -622,94 +550,25 @@ class OperatorCompiler(mappers.IdentityMapper):
self.assigned_names = set()
- # {{{ collecting various optemplate components ----------------------------
- def get_contained_fluxes(self, expr):
- """Recursively enumerate all flux expressions in the expression tree
- `expr`. The returned list consists of `ExecutionPlanner.FluxRecord`
- instances with fields `flux_expr` and `dependencies`.
- """
-
- contained_flux_ops = mappers.FluxCollector()(expr)
-
- from pytools import all
- assert all(isinstance(op, sym.WholeDomainFluxOperator)
- for op in contained_flux_ops), \
- "not all flux operators were of the expected type"
-
- return [self.FluxRecord(
- flux_expr=wdflux,
- dependencies=set(wdflux.interior_deps) | set(wdflux.boundary_deps),
- repr_op=wdflux.repr_op())
- for wdflux in contained_flux_ops]
+ # {{{ collect various optemplate components
def collect_diff_ops(self, expr):
- return mappers.BoundOperatorCollector(sym.ReferenceDiffOperatorBase)(expr)
-
- def collect_flux_exchange_ops(self, expr):
- return mappers.FluxExchangeCollector()(expr)
+ return mappers.BoundOperatorCollector(sym.RefDiffOperatorBase)(expr)
# }}}
- # {{{ top-level driver ----------------------------------------------------
- def __call__(self, expr, type_hints={}):
+ # {{{ top-level driver
+
+ def __call__(self, expr):
# Put the result expressions into variables as well.
expr = sym.cse(expr, "_result")
- from grudge.symbolic.mappers.type_inference import TypeInferrer
- self.typedict = TypeInferrer()(expr, type_hints)
-
- # {{{ flux batching
- # Fluxes can be evaluated faster in batches. Here, we find flux
- # batches that we can evaluate together.
-
- # For each FluxRecord, find the other fluxes its flux depends on.
- flux_queue = self.get_contained_fluxes(expr)
- for fr in flux_queue:
- fr.dependencies = set(sf.flux_expr
- for sf in self.get_contained_fluxes(fr.flux_expr)) \
- - set([fr.flux_expr])
-
- # Then figure out batches of fluxes to evaluate
- self.flux_batches = []
- admissible_deps = set()
- while flux_queue:
- present_batch = set()
- i = 0
- while i < len(flux_queue):
- fr = flux_queue[i]
- if fr.dependencies <= admissible_deps:
- present_batch.add(fr)
- flux_queue.pop(i)
- else:
- i += 1
-
- if present_batch:
- # bin batched operators by representative operator
- batches_by_repr_op = {}
- for fr in present_batch:
- batches_by_repr_op[fr.repr_op] = \
- batches_by_repr_op.get(fr.repr_op, set()) \
- | set([fr.flux_expr])
-
- for repr_op, batch in six.iteritems(batches_by_repr_op):
- self.flux_batches.append(
- self.FluxBatch(
- repr_op=repr_op,
- flux_exprs=list(batch)))
-
- admissible_deps |= set(fr.flux_expr for fr in present_batch)
- else:
- raise RuntimeError("cannot resolve flux evaluation order")
-
- # }}}
+ # from grudge.symbolic.mappers.type_inference import TypeInferrer
+ # self.typedict = TypeInferrer()(expr)
# Used for diff batching
-
self.diff_ops = self.collect_diff_ops(expr)
- # Flux exchange also works better when batched.
- self.flux_exchange_ops = self.collect_flux_exchange_ops(expr)
-
# Finally, walk the expression and build the code.
result = super(OperatorCompiler, self).__call__(expr)
@@ -719,7 +578,8 @@ class OperatorCompiler(mappers.IdentityMapper):
# }}}
- # {{{ variables and names -------------------------------------------------
+ # {{{ variables and names
+
def get_var_name(self, prefix=None):
def generate_suffixes():
yield ""
@@ -747,8 +607,7 @@ class OperatorCompiler(mappers.IdentityMapper):
self.assigned_names.add(name)
return name
- def assign_to_new_var(self, expr, priority=0, prefix=None,
- is_scalar_valued=False):
+ def assign_to_new_var(self, expr, priority=0, prefix=None):
from pymbolic.primitives import Variable, Subscript
# Observe that the only things that can be legally subscripted in
@@ -759,7 +618,7 @@ class OperatorCompiler(mappers.IdentityMapper):
new_name = self.get_var_name(prefix)
self.code.append(self.make_assign(
- new_name, expr, priority, is_scalar_valued))
+ new_name, expr, priority))
return Variable(new_name)
@@ -773,9 +632,6 @@ class OperatorCompiler(mappers.IdentityMapper):
except KeyError:
priority = getattr(expr, "priority", 0)
- from grudge.symbolic.mappers.type_inference import type_info
- is_scalar_valued = isinstance(self.typedict[expr], type_info.Scalar)
-
if isinstance(expr.child, sym.OperatorBinding):
# We need to catch operator bindings here and
# treat them specially. They get assigned to their
@@ -788,24 +644,14 @@ class OperatorCompiler(mappers.IdentityMapper):
rec_child = self.rec(expr.child)
cse_var = self.assign_to_new_var(rec_child,
- priority=priority, prefix=expr.prefix,
- is_scalar_valued=is_scalar_valued)
+ priority=priority, prefix=expr.prefix)
self.expr_to_var[expr.child] = cse_var
return cse_var
def map_operator_binding(self, expr, name_hint=None):
- from grudge.symbolic.operators import (
- ReferenceDiffOperatorBase,
- FluxOperatorBase)
-
- if isinstance(expr.op, ReferenceDiffOperatorBase):
+ if isinstance(expr.op, sym.RefDiffOperatorBase):
return self.map_ref_diff_op_binding(expr)
- elif isinstance(expr.op, FluxOperatorBase):
- raise RuntimeError("OperatorCompiler encountered a flux operator.\n\n"
- "We are expecting flux operators to be converted to custom "
- "flux assignment instructions, but the subclassed compiler "
- "does not seem to have done this.")
else:
# make sure operator assignments stand alone and don't get muddled
# up in vector math
@@ -859,15 +705,8 @@ class OperatorCompiler(mappers.IdentityMapper):
from pytools import single_valued
op_class = single_valued(type(d.op) for d in all_diffs)
- from grudge.symbolic.operators import \
- ReferenceQuadratureStiffnessTOperator
- if isinstance(op_class, ReferenceQuadratureStiffnessTOperator):
- assign_class = QuadratureDiffBatchAssign
- else:
- assign_class = DiffBatchAssign
-
self.code.append(
- assign_class(
+ DiffBatchAssign(
names=names,
op_class=op_class,
operators=[d.op for d in all_diffs],
@@ -880,84 +719,13 @@ class OperatorCompiler(mappers.IdentityMapper):
return self.expr_to_var[expr]
- def map_flux_exchange(self, expr):
- try:
- return self.expr_to_var[expr]
- except KeyError:
- from pytools.obj_array import obj_array_equal
- all_flux_xchgs = [fe
- for fe in self.flux_exchange_ops
- if obj_array_equal(fe.arg_fields, expr.arg_fields)]
-
- assert len(all_flux_xchgs) > 0
-
- names = [self.get_var_name() for d in all_flux_xchgs]
- self.code.append(
- FluxExchangeBatchAssign(
- names=names,
- indices_and_ranks=[
- (fe.index, fe.rank)
- for fe in all_flux_xchgs],
- arg_fields=[
- self.rec(arg_field)
- for arg_field in fe.arg_fields]))
-
- from pymbolic import var
- for n, d in zip(names, all_flux_xchgs):
- self.expr_to_var[d] = var(n)
-
- return self.expr_to_var[expr]
-
- def map_planned_flux(self, expr):
- try:
- return self.expr_to_var[expr]
- except KeyError:
- for fb in self.flux_batches:
- try:
- idx = fb.flux_exprs.index(expr)
- except ValueError:
- pass
- else:
- # found at idx
- mapped_fluxes = [
- self.internal_map_flux(f)
- for f in fb.flux_exprs]
-
- names = [self.get_var_name() for f in mapped_fluxes]
- self.code.append(
- self.make_flux_batch_assign(
- names, mapped_fluxes, fb.repr_op))
-
- from pymbolic import var
- for n, f in zip(names, fb.flux_exprs):
- self.expr_to_var[f] = var(n)
-
- return var(names[idx])
-
- raise RuntimeError("flux '%s' not in any flux batch" % expr)
-
# }}}
# {{{ instruction producers
- def make_assign(self, name, expr, priority, is_scalar_valued=False):
+ def make_assign(self, name, expr, priority):
return Assign(names=[name], exprs=[expr],
- priority=priority,
- is_scalar_valued=is_scalar_valued)
-
- def make_flux_batch_assign(self, names, expressions, repr_op):
- return FluxBatchAssign(names=names, expressions=expressions, repr_op=repr_op)
-
- # from CUDA backend:
- # def make_flux_batch_assign(self, names, expressions, repr_op):
- # from pytools import single_valued
- # quadrature_tag = single_valued(
- # wdflux.quadrature_tag
- # for wdflux in expressions)
-
- # return CUDAFluxBatchAssign(
- # names=names, expressions=expressions, repr_op=repr_op,
- # quadrature_tag=quadrature_tag)
+ priority=priority)
# }}}
@@ -1164,22 +932,6 @@ class OperatorCompiler(mappers.IdentityMapper):
# }}}
- def internal_map_flux(self, wdflux):
- return sym.WholeDomainFluxOperator(
- wdflux.is_lift,
- [wdflux.InteriorInfo(
- flux_expr=ii.flux_expr,
- field_expr=self.rec(ii.field_expr))
- for ii in wdflux.interiors],
- [wdflux.BoundaryInfo(
- flux_expr=bi.flux_expr,
- bpair=self.rec(bi.bpair))
- for bi in wdflux.boundaries],
- wdflux.quadrature_tag)
-
- def map_whole_domain_flux(self, wdflux):
- return self.map_planned_flux(wdflux)
-
def finalize_multi_assign(self, names, exprs, do_not_return, priority):
return VectorExprAssign(names=names, exprs=exprs,
do_not_return=do_not_return,
diff --git a/grudge/symbolic/flux/__init__.py b/grudge/symbolic/flux/__init__.py
deleted file mode 100644
index e69de29b..00000000
diff --git a/grudge/symbolic/flux/mappers.py b/grudge/symbolic/flux/mappers.py
deleted file mode 100644
index 1533274a..00000000
--- a/grudge/symbolic/flux/mappers.py
+++ /dev/null
@@ -1,188 +0,0 @@
-"""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.py b/grudge/symbolic/flux/primitives.py
deleted file mode 100644
index e2629fbe..00000000
--- a/grudge/symbolic/flux/primitives.py
+++ /dev/null
@@ -1,256 +0,0 @@
-"""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 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
deleted file mode 100644
index baaee7d5..00000000
--- a/grudge/symbolic/flux/tools.py
+++ /dev/null
@@ -1,85 +0,0 @@
-"""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
index 815e6244..06518a3b 100644
--- a/grudge/symbolic/mappers/__init__.py
+++ b/grudge/symbolic/mappers/__init__.py
@@ -38,7 +38,6 @@ import pymbolic.mapper.flop_counter
from pymbolic.mapper import CSECachingMapperMixin
from grudge import sym
-from grudge import sym_flux
import grudge.symbolic.operators as op
@@ -128,29 +127,34 @@ class OperatorReducerMixin(LocalOpReducerMixin, FluxOpReducerMixin):
def _map_op_base(self, expr, *args, **kwargs):
return self.map_operator(expr, *args, **kwargs)
+ map_elementwise_linear = _map_op_base
+
+ map_interpolation = _map_op_base
+
map_nodal_sum = _map_op_base
- map_nodal_max = _map_op_base
map_nodal_min = _map_op_base
+ map_nodal_max = _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
+ map_stiffness = _map_op_base
+ map_diff = _map_op_base
+ map_stiffness_t = _map_op_base
+
+ map_ref_diff = _map_op_base
+ map_ref_stiffness_t = _map_op_base
+
+ map_mass = _map_op_base
+ map_inverse_mass = _map_op_base
+ map_ref_mass = _map_op_base
+ map_ref_inverse_mass = _map_op_base
+
+ map_opposite_interior_face_swap = _map_op_base
+ map_face_mass_operator = _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):
@@ -158,19 +162,11 @@ class IdentityMapperMixin(LocalOpReducerMixin, FluxOpReducerMixin):
"IdentityMapper instances cannot be combined with " \
"the BoundOpMapperMixin"
- return expr.__class__(
+ return type(expr)(
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)
+ # {{{ operators
def map_elementwise_linear(self, expr, *args, **kwargs):
assert not isinstance(self, BoundOpMapperMixin), \
@@ -180,30 +176,41 @@ class IdentityMapperMixin(LocalOpReducerMixin, FluxOpReducerMixin):
# 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_interpolation = map_elementwise_linear
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
+ map_stiffness = map_elementwise_linear
+ map_diff = map_elementwise_linear
+ map_stiffness_t = map_elementwise_linear
+
+ map_ref_diff = map_elementwise_linear
+ map_ref_stiffness_t = map_elementwise_linear
+
+ map_mass = map_elementwise_linear
+ map_inverse_mass = map_elementwise_linear
+ map_ref_mass = map_elementwise_linear
+ map_ref_inverse_mass = map_elementwise_linear
+
+ map_opposite_interior_face_swap = map_elementwise_linear
+ map_face_mass_operator = map_elementwise_linear
+
+ # }}}
+
+ # {{{ primitives
+
+ def map_grudge_variable(self, expr, *args, **kwargs):
+ # it's a leaf--no changing children
+ return expr
+
+ map_c_function = map_grudge_variable
+
+ map_ones = map_grudge_variable
+ map_node_coordinate_component = map_grudge_variable
+
+ # }}}
class BoundOpMapperMixin(object):
@@ -247,14 +254,13 @@ class DependencyMapper(
def map_operator(self, expr):
return set()
- def map_scalar_parameter(self, expr):
+ def map_grudge_variable(self, expr):
return set([expr])
def _map_leaf(self, expr):
return set()
map_ones = _map_leaf
-
map_node_coordinate_component = _map_leaf
@@ -377,12 +383,13 @@ class OperatorSpecializer(CSECachingMapperMixin, IdentityMapper):
# operators.
# {{{ elementwise operators
+
if isinstance(expr.op, op.MassOperator) and has_quad_operand:
return op.QuadratureMassOperator(
field_repr_tag.quadrature_tag)(self.rec(expr.field))
- elif isinstance(expr.op, op.ReferenceMassOperator) and has_quad_operand:
- return op.ReferenceQuadratureMassOperator(
+ elif isinstance(expr.op, op.RefMassOperator) and has_quad_operand:
+ return op.RefQuadratureMassOperator(
field_repr_tag.quadrature_tag)(self.rec(expr.field))
elif (isinstance(expr.op, op.StiffnessTOperator) and has_quad_operand):
@@ -390,9 +397,9 @@ class OperatorSpecializer(CSECachingMapperMixin, IdentityMapper):
expr.op.xyz_axis, field_repr_tag.quadrature_tag)(
self.rec(expr.field))
- elif (isinstance(expr.op, op.ReferenceStiffnessTOperator)
+ elif (isinstance(expr.op, op.RefStiffnessTOperator)
and has_quad_operand):
- return op.ReferenceQuadratureStiffnessTOperator(
+ return op.RefQuadratureStiffnessTOperator(
expr.op.xyz_axis, field_repr_tag.quadrature_tag)(
self.rec(expr.field))
@@ -489,20 +496,12 @@ class GlobalToReferenceMapper(CSECachingMapperMixin, IdentityMapper):
# if we encounter non-quadrature operators here, we know they
# must be nodal.
- jac_notag = sym.area_element(self.ambient_dim, self.dim,
- where=expr.op.where, quadrature_tag=None)
-
- def rewrite_derivative(ref_class, field, where,
- input_quadrature_tag=None, with_jacobian=True):
- jac_tag = sym.area_element(self.ambient_dim, self.dim,
- where=expr.op.where, quadrature_tag=input_quadrature_tag)
-
- if input_quadrature_tag is not None:
- diff_kwargs = dict(input_quadrature_tag=input_quadrature_tag)
- else:
- diff_kwargs = {}
+ jac_in = sym.area_element(self.ambient_dim, self.dim, dd=expr.op.dd_in)
+ jac_noquad = sym.area_element(self.ambient_dim, self.dim,
+ dd=expr.op.dd_in.with_qtag(sym.QTAG_NONE))
- diff_kwargs["where"] = where
+ def rewrite_derivative(ref_class, field, dd_in, with_jacobian=True):
+ jac_tag = sym.area_element(self.ambient_dim, self.dim, dd=dd_in)
rec_field = self.rec(field)
if with_jacobian:
@@ -511,45 +510,31 @@ class GlobalToReferenceMapper(CSECachingMapperMixin, IdentityMapper):
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)
+ ref_class(rst_axis, dd_in=dd_in)(rec_field)
for rst_axis in range(self.dim))
if isinstance(expr.op, op.MassOperator):
- return op.ReferenceMassOperator()(
- 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, expr.op.where)(
- jac_tag * self.rec(expr.field))
+ return op.RefMassOperator(expr.op.dd_in, expr.op.dd_out)(
+ jac_in * self.rec(expr.field))
elif isinstance(expr.op, op.InverseMassOperator):
- return op.ReferenceInverseMassOperator(expr.op.where)(
- 1/jac_notag * self.rec(expr.field))
+ return op.RefInverseMassOperator(expr.op.dd_in, expr.op.dd_out)(
+ 1/jac_in * self.rec(expr.field))
elif isinstance(expr.op, op.StiffnessOperator):
- return op.ReferenceMassOperator(expr.op.where)(jac_notag *
+ return op.RefMassOperator()(jac_noquad *
self.rec(
- op.DifferentiationOperator(expr.op.xyz_axis)(expr.field)))
+ op.DiffOperator(expr.op.xyz_axis)(expr.field)))
- elif isinstance(expr.op, op.DifferentiationOperator):
+ elif isinstance(expr.op, op.DiffOperator):
return rewrite_derivative(
- op.ReferenceDifferentiationOperator,
- expr.field, where=expr.op.where, with_jacobian=False)
+ op.RefDiffOperator,
+ expr.field, dd_in=expr.op.dd_in, with_jacobian=False)
elif isinstance(expr.op, op.StiffnessTOperator):
return rewrite_derivative(
- op.ReferenceStiffnessTOperator,
- expr.field, where=expr.op.where)
-
- elif isinstance(expr.op, op.QuadratureStiffnessTOperator):
- return rewrite_derivative(
- op.ReferenceQuadratureStiffnessTOperator,
- expr.field, input_quadrature_tag=expr.op.input_quadrature_tag,
- where=expr.op.where)
+ op.RefStiffnessTOperator,
+ expr.field, dd_in=expr.op.dd_in)
elif isinstance(expr.op, op.MInvSTOperator):
return self.rec(
@@ -565,72 +550,67 @@ class GlobalToReferenceMapper(CSECachingMapperMixin, IdentityMapper):
# {{{ stringification ---------------------------------------------------------
class StringifyMapper(pymbolic.mapper.stringifier.StringifyMapper):
- def _format_btag(self, tag):
- if isinstance(tag, type):
- return tag.__name__
+ def _format_dd(self, dd):
+ def fmt(s):
+ if isinstance(s, type):
+ return s.__name__
+ else:
+ return repr(s)
+
+ from meshmode.discretization.connection import (
+ FRESTR_ALL_FACES, FRESTR_INTERIOR_FACES)
+ if dd.domain_tag is None:
+ result = "?"
+ elif dd.domain_tag is sym.DTAG_VOLUME_ALL:
+ result = "vol"
+ elif dd.domain_tag is sym.DTAG_SCALAR:
+ result = "scalar"
+ elif dd.domain_tag is FRESTR_ALL_FACES:
+ result = "all_faces"
+ elif dd.domain_tag is FRESTR_INTERIOR_FACES:
+ result = "int_faces"
else:
- return repr(tag)
+ result = fmt(dd.domain_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
+ if dd.quadrature_tag is None:
+ pass
+ elif dd.quadrature_tag is sym.QTAG_NONE:
+ result += "q"
+ else:
+ result += "Q"+fmt(dd.quadrature_tag)
return result
- 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.symbolic.flux.mappers 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))
+ def _format_op_dd(self, op):
+ return "[%s->%s]" % (self._format_dd(op.dd_in), self._format_dd(op.dd_out))
# {{{ nodal ops
def map_nodal_sum(self, expr, enclosing_prec):
- return "NodalSum"
+ return "NodalSum" + self._format_op_dd(expr)
def map_nodal_max(self, expr, enclosing_prec):
- return "NodalMax"
+ return "NodalMax" + self._format_op_dd(expr)
def map_nodal_min(self, expr, enclosing_prec):
- return "NodalMin"
+ return "NodalMin" + self._format_op_dd(expr)
# }}}
# {{{ global differentiation
def map_diff(self, expr, enclosing_prec):
- return "Diffx%d%s" % (expr.xyz_axis, self._format_op_props(expr))
+ return "Diffx%d%s" % (expr.xyz_axis, self._format_op_dd(expr))
def map_minv_st(self, expr, enclosing_prec):
- return "MInvSTx%d%s" % (expr.xyz_axis, self._format_op_props(expr))
+ return "MInvSTx%d%s" % (expr.xyz_axis, self._format_op_dd(expr))
def map_stiffness(self, expr, enclosing_prec):
- return "Stiffx%d%s" % (expr.xyz_axis, self._format_op_props(expr))
+ return "Stiffx%d%s" % (expr.xyz_axis, self._format_op_dd(expr))
def map_stiffness_t(self, expr, enclosing_prec):
- return "StiffTx%d%s" % (expr.xyz_axis, self._format_op_props(expr))
+ return "StiffTx%d%s" % (expr.xyz_axis, self._format_op_dd(expr))
- def map_quad_stiffness_t(self, expr, enclosing_prec):
- return "StiffTx%d%s" % (
- expr.input_quadrature_tag, expr.xyz_axis,
- self._format_op_props(expr))
# }}}
# {{{ global mass
@@ -641,107 +621,49 @@ class StringifyMapper(pymbolic.mapper.stringifier.StringifyMapper):
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%s" % (expr.rst_axis, self._format_op_props(expr))
+ return "Diffr%d%s" % (expr.rst_axis, self._format_op_dd(expr))
def map_ref_stiffness_t(self, expr, enclosing_prec):
- return "StiffTr%d%s" % (expr.rst_axis, self._format_op_props(expr))
+ return "StiffTr%d%s" % (expr.rst_axis, self._format_op_dd(expr))
- def map_ref_quad_stiffness_t(self, expr, enclosing_prec):
- 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" + self._format_op_props(expr)
+ return "RefM" + self._format_op_dd(expr)
def map_ref_inverse_mass(self, expr, enclosing_prec):
- return "RefInvM" + self._format_op_props(expr)
-
- def map_ref_quad_mass(self, expr, enclosing_prec):
- return "RefM" + self._format_op_props(expr)
+ return "RefInvM" + self._format_op_dd(expr)
# }}}
def map_elementwise_linear(self, expr, enclosing_prec):
return "ElWLin:%s%s" % (
expr.__class__.__name__,
- self._format_op_props(expr))
+ self._format_op_dd(expr))
# {{{ 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_face_mass_operator(self, expr, enclosing_prec):
+ return "FaceM" + self._format_op_dd(expr)
- 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_opposite_interior_face_swap(self, expr, enclosing_prec):
+ return "OppSwap" + self._format_op_dd(expr)
- 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" + self._format_op_props(expr)
-
- def map_boundarize(self, expr, enclosing_prec):
- 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
- 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):
return "Ones" + self._format_op_props(expr)
# {{{ 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))
+ return "x[%d]@%s" % (expr.axis, self._format_dd(expr.dd))
# }}}
@@ -754,100 +676,17 @@ class StringifyMapper(pymbolic.mapper.stringifier.StringifyMapper):
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 "Upsamp" + self._format_op_props(expr)
+ def map_grudge_variable(self, expr, enclosing_prec):
+ return "%s:%s" % (expr.name, self._format_dd(expr.dd))
- def map_quad_int_faces_grid_upsampler(self, expr, enclosing_prec):
- return "ToIntFaceQ[%s]" % expr.quadrature_tag
+ def map_interpolation(self, expr, enclosing_prec):
+ return "Interp" + self._format_op_dd(expr)
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.symbolic.flux.mappers 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, self._format_btag(expr.tag))
-
- def map_operator_binding(self, expr, enclosing_prec):
- if isinstance(expr.op, op.RestrictToBoundary):
- from pymbolic.mapper.stringifier import PREC_CALL, PREC_SUM
- return self.parenthesize_if_needed(
- "%s@%s" % (
- self.rec(expr.field, PREC_CALL),
- self._format_btag(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))
+ pass
class NoCSEStringifyMapper(StringifyMapper):
@@ -974,115 +813,7 @@ class CommutativeConstantFoldingMapper(
if is_zero(field):
return 0
- if isinstance(expr.op, op.FluxOperatorBase):
- if isinstance(field, sym.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, make_obj_array
- if not is_obj_array(vol_field):
- vol_field = [vol_field]
-
- subst_map = {}
-
- from grudge.tools import is_zero
-
- # 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[sym_flux.FieldComponent(i, is_interior=True)] = 0
- subst_map[sym_flux.FieldComponent(i, is_interior=False)] = 0
- else:
- new_vol_field.append(flux_arg)
- subst_map[sym_flux.FieldComponent(i, is_interior=True)] = \
- sym_flux.FieldComponent(new_idx, is_interior=True)
- subst_map[sym_flux.FieldComponent(i, is_interior=False)] = \
- sym_flux.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.symbolic.flux.mappers 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, make_obj_array
- if not is_obj_array(vol_field):
- vol_field = [vol_field]
- if not is_obj_array(bdry_field):
- bdry_field = [bdry_field]
-
- subst_map = {}
-
- # process volume field
- from grudge.tools import is_zero
- new_vol_field = []
- new_idx = 0
- for i, flux_arg in enumerate(vol_field):
- fc = sym_flux.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] = sym_flux.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 = sym_flux.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] = sym_flux.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.symbolic.flux.mappers import FluxSubstitutionMapper
- new_flux = FluxSubstitutionMapper(sub_flux)(op.flux)
-
- if is_zero(new_flux):
- return 0
- else:
- from grudge.symbolic.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))
-
- # }}}
+ return expr.op(field)
class EmptyFluxKiller(CSECachingMapperMixin, IdentityMapper):
@@ -1094,8 +825,13 @@ class EmptyFluxKiller(CSECachingMapperMixin, IdentityMapper):
IdentityMapper.map_common_subexpression
def map_operator_binding(self, expr):
- if (isinstance(expr.op, op.BoundaryFluxOperatorBase) and
- len(self.mesh.tag_to_boundary.get(expr.op.boundary_tag, [])) == 0):
+ from meshmode.mesh import is_boundary_tag_empty
+ if (isinstance(expr.op, sym.InterpolationOperator)
+ and expr.op.dd_out.is_boundary()
+ and expr.op.dd_out.domain_tag not in [
+ sym.FRESTR_ALL_FACES, sym.FRESTR_INTERIOR_FACES]
+ and is_boundary_tag_empty(self.mesh,
+ expr.op.dd_out.domain_tag)):
return 0
else:
return IdentityMapper.map_operator_binding(self, expr)
@@ -1139,7 +875,7 @@ class _InnerDerivativeJoiner(pymbolic.mapper.RecursiveMapper):
else:
return self.rec(expr, derivatives)
- for operator, operands in six.iteritem(sub_derivatives):
+ for operator, operands in six.iteritems(sub_derivatives):
for operand in operands:
derivatives.setdefault(operator, []).append(
factor*operand)
@@ -1283,9 +1019,6 @@ class InverseMassContractor(CSECachingMapperMixin, IdentityMapper):
map_common_subexpression_uncached = \
IdentityMapper.map_common_subexpression
- def map_boundary_pair(self, bp):
- return sym.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
@@ -1333,62 +1066,23 @@ class ErrorChecker(CSECachingMapperMixin, IdentityMapper):
# }}}
-# {{{ collectors for various optemplate components
+# {{{ collectors for various symbolic operator components
class CollectorMixin(OperatorReducerMixin, LocalOpReducerMixin, FluxOpReducerMixin):
def combine(self, values):
from pytools import flatten
return set(flatten(values))
- def map_constant(self, bpair):
+ def map_constant(self, expr, *args, **kwargs):
return set()
- map_operator = map_constant
- map_variable = map_constant
- map_ones = map_constant
- map_node_coordinate_component = 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):
- if isinstance(expr.op, 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
+ map_grudge_variable = map_constant
+ map_c_function = map_grudge_variable
+ map_ones = map_grudge_variable
+ map_node_coordinate_component = map_grudge_variable
-class BoundaryTagCollector(CollectorMixin, CombineMapper):
- def map_boundary_pair(self, bpair):
- return set([bpair.tag])
+ map_operator = map_grudge_variable
# I'm not sure this works still.
@@ -1425,113 +1119,7 @@ class FluxExchangeCollector(CSECachingMapperMixin, CollectorMixin, CombineMapper
# {{{ evaluation
class Evaluator(pymbolic.mapper.evaluator.EvaluationMapper):
- def map_boundary_pair(self, bp):
- return sym.BoundaryPair(self.rec(bp.field), self.rec(bp.bfield), bp.tag)
-# }}}
-
-
-# {{{ boundary combiner
-
-class BoundaryCombiner(CSECachingMapperMixin, IdentityMapper):
- """Combines inner fluxes and boundary fluxes into a
- single, whole-domain operator of type
- :class:`grudge.symbolic.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.symbolic.primitives import OperatorBinding, BoundaryPair
-
- 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:
- if (isinstance(ch, OperatorBinding)
- and isinstance(ch.op, 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
-
- if isinstance(ch.op, 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)
- from meshmode.mesh import is_boundary_tag_empty
- if not is_boundary_tag_empty(self.mesh, bpair.tag):
- boundaries.append(op.WholeDomainFluxOperator.BoundaryInfo(
- flux_expr=ch.op.flux,
- bpair=bpair))
- else:
- interiors.append(op.WholeDomainFluxOperator.InteriorInfo(
- flux_expr=ch.op.flux,
- field_expr=self.rec(ch.field)))
- else:
- rest.append(ch)
-
- if interiors or boundaries:
- wdf = op.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):
- if isinstance(expr.op, 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)
+ pass
# }}}
diff --git a/grudge/symbolic/operators.py b/grudge/symbolic/operators.py
index de96dba9..631fac36 100644
--- a/grudge/symbolic/operators.py
+++ b/grudge/symbolic/operators.py
@@ -27,23 +27,36 @@ THE SOFTWARE.
from six.moves import intern
import numpy as np
-import numpy.linalg as la
+import numpy.linalg as la # noqa
import pymbolic.primitives
-from pytools import Record, memoize_method
+
+
+def _sym():
+ # A hack to make referring to grudge.sym possible without
+ # circular imports and tons of typing.
+
+ from grudge import sym
+ return sym
# {{{ base classes
-class Operator(pymbolic.primitives.Leaf):
+class Operator(pymbolic.primitives.Expression):
"""
- .. attribute:: where
+ .. attribute:: dd_in
+
+ an instance of :class:`grudge.sym.DOFDesc` describing the
+ input discretization.
+
+ .. attribute:: dd_out
- *None* for the default volume discretization or a boundary
- tag for an operation on the denoted part of the boundary.
+ an instance of :class:`grudge.sym.DOFDesc` describing the
+ output discretization.
"""
- def __init__(self, where=None):
- self.where = where
+ def __init__(self, dd_in, dd_out):
+ self.dd_in = _sym().as_dofdesc(dd_in)
+ self.dd_out = _sym().as_dofdesc(dd_out)
def stringifier(self):
from grudge.symbolic.mappers import StringifyMapper
@@ -62,23 +75,47 @@ class Operator(pymbolic.primitives.Leaf):
return with_object_array_or_scalar(bind_one, expr)
- def get_hash(self):
- return hash((self.__class__,) + (self.__getinitargs__()))
-
- def is_equal(self, other):
- return self.__class__ == other.__class__ and \
- self.__getinitargs__() == other.__getinitargs__()
+ def with_dd(self, dd_in=None, dd_out=None):
+ """Return a copy of *self*, modified to the given DOF descriptors.
+ """
+ return type(self)(
+ *self.__getinitargs__()[:-2],
+ dd_in=dd_in or self.dd_in,
+ dd_out=dd_out or self.dd_out)
def __getinitargs__(self):
- return (self.where,)
+ return (self.dd_in, self.dd_out,)
# }}}
-# {{{ sum, integral, max
+class ElementwiseLinearOperator(Operator):
+ def matrix(self, element_group):
+ raise NotImplementedError
+
+ mapper_method = intern("map_elementwise_linear")
+
+
+class InterpolationOperator(Operator):
+ mapper_method = intern("map_interpolation")
+
+interp = InterpolationOperator
+
+
+class ElementwiseMaxOperator(Operator):
+ mapper_method = intern("map_elementwise_max")
+
+
+# {{{ nodal reduction: sum, integral, max
class NodalReductionOperator(Operator):
- pass
+ def __init__(self, dd_in, dd_out=None):
+ if dd_out is None:
+ dd_out = _sym().DD_SCALAR
+
+ assert dd_out.is_scalar()
+
+ super(NodalReductionOperator, self).__init__(dd_out=dd_out, dd_in=dd_in)
class NodalSum(NodalReductionOperator):
@@ -100,16 +137,18 @@ class NodalMin(NodalReductionOperator):
# {{{ global differentiation
class DiffOperatorBase(Operator):
- def __init__(self, xyz_axis):
- Operator.__init__(self)
+ def __init__(self, xyz_axis, dd_in=None, dd_out=None):
+ if dd_in is None:
+ dd_in = _sym().DD_VOLUME
+ if dd_out is None:
+ dd_out = dd_in.with_qtag(_sym().QTAG_NONE)
+
+ super(DiffOperatorBase, self).__init__(dd_in, dd_out)
self.xyz_axis = xyz_axis
def __getinitargs__(self):
- return (self.xyz_axis,)
-
- def preimage_ranges(self, eg):
- return eg.ranges
+ return (self.xyz_axis, self.dd_in, self.dd_out)
def equal_except_for_axis(self, other):
return (type(self) == type(other)
@@ -129,7 +168,7 @@ class StiffnessOperator(StrongFormDiffOperatorBase):
mapper_method = intern("map_stiffness")
-class DifferentiationOperator(StrongFormDiffOperatorBase):
+class DiffOperator(StrongFormDiffOperatorBase):
mapper_method = intern("map_diff")
@@ -140,46 +179,24 @@ class StiffnessTOperator(WeakFormDiffOperatorBase):
class MInvSTOperator(WeakFormDiffOperatorBase):
mapper_method = intern("map_minv_st")
-
-class QuadratureStiffnessTOperator(DiffOperatorBase):
- """
- .. note::
-
- This operator is purely for internal use. It is inserted
- by :class:`grudge.symbolic.mappers.OperatorSpecializer`
- when a :class:`StiffnessTOperator` is applied to a quadrature
- field, and then eliminated by
- :class:`grudge.symbolic.mappers.GlobalToReferenceMapper`
- in favor of operators on the reference element.
- """
-
- 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.input_quadrature_tag, self.where)
-
- mapper_method = intern("map_quad_stiffness_t")
-
-
-def DiffOperatorVector(els):
- from grudge.tools import join_fields
- return join_fields(*els)
-
# }}}
# {{{ reference-element differentiation
-class ReferenceDiffOperatorBase(Operator):
- def __init__(self, rst_axis, where=None):
- super(ReferenceDiffOperatorBase, self).__init__(where)
+class RefDiffOperatorBase(ElementwiseLinearOperator):
+ def __init__(self, rst_axis, dd_in=None, dd_out=None):
+ if dd_in is None:
+ dd_in = _sym().DD_VOLUME
+ if dd_out is None:
+ dd_out = dd_in.with_qtag(_sym().QTAG_NONE)
+
+ super(RefDiffOperatorBase, self).__init__(dd_in, dd_out)
self.rst_axis = rst_axis
def __getinitargs__(self):
- return (self.rst_axis, self.where)
+ return (self.rst_axis, self.dd_in, self.dd_out)
def equal_except_for_axis(self, other):
return (type(self) == type(other)
@@ -187,93 +204,22 @@ class ReferenceDiffOperatorBase(Operator):
and self.__getinitargs__()[1:] == other.__getinitargs__()[1:])
-class ReferenceDifferentiationOperator(ReferenceDiffOperatorBase):
+class RefDiffOperator(RefDiffOperatorBase):
mapper_method = intern("map_ref_diff")
-class ReferenceStiffnessTOperator(ReferenceDiffOperatorBase):
+class RefStiffnessTOperator(RefDiffOperatorBase):
mapper_method = intern("map_ref_stiffness_t")
-
-class ReferenceQuadratureStiffnessTOperator(ReferenceDiffOperatorBase):
- """
- .. note::
-
- This operator is purely for internal use. It is inserted
- by :class:`grudge.symbolic.mappers.OperatorSpecializer`
- when a :class:`StiffnessTOperator` is applied to a quadrature field.
- """
-
- 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.input_quadrature_tag,
- self.where)
-
- mapper_method = intern("map_ref_quad_stiffness_t")
-
-# }}}
-
# }}}
-
-# {{{ elementwise operators
-
-class ElementwiseLinearOperator(Operator):
- def matrix(self, element_group):
- raise NotImplementedError
-
- mapper_method = intern("map_elementwise_linear")
-
-
-class ElementwiseMaxOperator(Operator):
- mapper_method = intern("map_elementwise_max")
-
-
-# {{{ quadrature upsamplers
-
-class QuadratureGridUpsampler(Operator):
- """In a user-specified optemplate, this operator can be used to interpolate
- volume and boundary data to their corresponding quadrature grids.
-
- In pre-processing, the boundary quad interpolation is specialized to
- a separate operator, :class:`QuadratureBoundaryGridUpsampler`.
- """
- def __init__(self, quadrature_tag, where=None):
- self.quadrature_tag = quadrature_tag
- self.where = where
-
- def __getinitargs__(self):
- return (self.quadrature_tag,)
-
- mapper_method = intern("map_quad_grid_upsampler")
-
-
-class QuadratureInteriorFacesGridUpsampler(Operator):
- """Interpolates nodal volume data to interior face data on a quadrature
- grid.
-
- Note that the "interior faces" grid includes faces lying opposite to the
- boundary.
- """
- def __init__(self, quadrature_tag):
- self.quadrature_tag = quadrature_tag
-
- def __getinitargs__(self):
- return (self.quadrature_tag,)
-
- mapper_method = intern("map_quad_int_faces_grid_upsampler")
-
# }}}
# {{{ various elementwise linear operators
class FilterOperator(ElementwiseLinearOperator):
- def __init__(self, mode_response_func):
+ def __init__(self, mode_response_func, dd_in=None, dd_out=None):
"""
:param mode_response_func: A function mapping
``(mode_tuple, local_discretization)`` to a float indicating the
@@ -281,12 +227,27 @@ class FilterOperator(ElementwiseLinearOperator):
(For example an instance of
:class:`ExponentialFilterResponseFunction`.
"""
+ if dd_in is None:
+ dd_in = _sym().DD_VOLUME
+ if dd_out is None:
+ dd_out = dd_in
+
+ if dd_in.uses_quadrature():
+ raise ValueError("dd_in may not use quadrature")
+ if dd_in != dd_out:
+ raise ValueError("dd_in and dd_out must be identical")
+
+ super(FilterOperator, self).__init__(dd_in, dd_out)
+
self.mode_response_func = mode_response_func
def __getinitargs__(self):
- return (self.mode_response_func,)
+ return (self.mode_response_func, self.dd_in, self.dd_out)
def matrix(self, eg):
+ # FIXME
+ raise NotImplementedError()
+
ldis = eg.local_discretization
filter_coeffs = [self.mode_response_func(mid, ldis)
@@ -303,19 +264,27 @@ class FilterOperator(ElementwiseLinearOperator):
return mat
-class OnesOperator(ElementwiseLinearOperator):
- def matrix(self, eg):
- ldis = eg.local_discretization
+class AveragingOperator(ElementwiseLinearOperator):
+ def __init__(self, dd_in=None, dd_out=None):
+ if dd_in is None:
+ dd_in = _sym().DD_VOLUME
+ if dd_out is None:
+ dd_out = dd_in
- node_count = ldis.node_count()
- return np.ones((node_count, node_count), dtype=np.float64)
+ if dd_in.uses_quadrature():
+ raise ValueError("dd_in may not use quadrature")
+ if dd_in != dd_out:
+ raise ValueError("dd_in and dd_out must be identical")
+ super(FilterOperator, self).__init__(dd_in, dd_out)
-class AveragingOperator(ElementwiseLinearOperator):
def matrix(self, eg):
# average matrix, so that AVE*fields = cellaverage(fields)
# see Hesthaven and Warburton page 227
+ # FIXME
+ raise NotImplementedError()
+
mmat = eg.local_discretization.mass_matrix()
standard_el_vol = np.sum(np.dot(mmat, np.ones(mmat.shape[0])))
avg_mat_row = np.sum(mmat, 0)/standard_el_vol
@@ -326,17 +295,13 @@ class AveragingOperator(ElementwiseLinearOperator):
class InverseVandermondeOperator(ElementwiseLinearOperator):
- def matrix(self, eg):
- return np.asarray(
- la.inv(eg.local_discretization.vandermonde()),
- order="C")
+ def matrix(self, element_group):
+ raise NotImplementedError() # FIXME
class VandermondeOperator(ElementwiseLinearOperator):
- def matrix(self, eg):
- return np.asarray(
- eg.local_discretization.vandermonde(),
- order="C")
+ def matrix(self, element_group):
+ raise NotImplementedError() # FIXME
# }}}
@@ -346,7 +311,19 @@ class VandermondeOperator(ElementwiseLinearOperator):
# {{{ mass operators
class MassOperatorBase(Operator):
- pass
+ def __init__(self, dd_in=None, dd_out=None):
+ if dd_in is None:
+ dd_in = _sym().DD_VOLUME
+ if dd_out is None:
+ dd_out = dd_in
+
+ if not dd_in.is_volume():
+ raise ValueError("mass operators only work on volume "
+ "discretization")
+ if dd_out != dd_in:
+ raise ValueError("dd_out and dd_in must be identical")
+
+ super(MassOperatorBase, self).__init__(dd_in, dd_out)
class MassOperator(MassOperatorBase):
@@ -357,51 +334,11 @@ class InverseMassOperator(MassOperatorBase):
mapper_method = intern("map_inverse_mass")
-class QuadratureMassOperator(Operator):
- """
- .. note::
-
- This operator is purely for internal use. It is inserted
- by :class:`grudge.symbolic.mappers.OperatorSpecializer`
- when a :class:`StiffnessTOperator` is applied to a quadrature
- field, and then eliminated by
- :class:`grudge.symbolic.mappers.GlobalToReferenceMapper`
- in favor of operators on the reference element.
- """
-
- def __init__(self, quadrature_tag):
- self.quadrature_tag = quadrature_tag
-
- def __getinitargs__(self):
- return (self.quadrature_tag,)
-
- mapper_method = intern("map_quad_mass")
-
-
-class ReferenceQuadratureMassOperator(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, where=None):
- super(Operator, self).__init__(self.where)
- self.quadrature_tag = quadrature_tag
-
- def __getinitargs__(self):
- return (self.quadrature_tag, self.where)
-
- mapper_method = intern("map_ref_quad_mass")
-
-
-class ReferenceMassOperatorBase(MassOperatorBase):
+class RefMassOperatorBase(ElementwiseLinearOperator):
pass
-class ReferenceMassOperator(ReferenceMassOperatorBase):
+class RefMassOperator(RefMassOperatorBase):
@staticmethod
def matrix(element_group):
import modepy as mp
@@ -412,7 +349,7 @@ class ReferenceMassOperator(ReferenceMassOperatorBase):
mapper_method = intern("map_ref_mass")
-class ReferenceInverseMassOperator(ReferenceMassOperatorBase):
+class RefInverseMassOperator(RefMassOperatorBase):
@staticmethod
def matrix(element_group):
import modepy as mp
@@ -427,319 +364,53 @@ class ReferenceInverseMassOperator(ReferenceMassOperatorBase):
# {{{ boundary-related operators
-class RestrictToBoundary(Operator):
- def __init__(self, tag):
- self.tag = tag
-
- def __getinitargs__(self):
- return (self.tag,)
-
- mapper_method = intern("map_boundarize")
-
-
-class FluxExchangeOperator(pymbolic.primitives.AlgebraicLeaf):
- """An operator that results in the sending and receiving of
- boundary information for its argument fields.
- """
-
- def __init__(self, idx, rank, arg_fields):
- self.index = idx
- self.rank = rank
- self.arg_fields = arg_fields
-
- # only tuples are hashable
- if not isinstance(arg_fields, tuple):
- raise TypeError("FluxExchangeOperator: arg_fields must be a tuple")
-
- def __getinitargs__(self):
- return (self.index, self.rank, self.arg_fields)
-
- def get_hash(self):
- return hash((self.__class__, self.index, self.rank, self.arg_fields))
-
- mapper_method = intern("map_flux_exchange")
+class OppositeInteriorFaceSwap(Operator):
+ def __init__(self, dd_in=None, dd_out=None):
+ sym = _sym()
- def is_equal(self, other):
- return self.__class__ == other.__class__ and \
- self.__getinitargs__() == other.__getinitargs__()
+ if dd_in is None:
+ dd_in = sym.DOFDesc(sym.FRESTR_INTERIOR_FACES, None)
+ if dd_out is None:
+ dd_out = dd_in
-# }}}
-
-
-# {{{ flux-like operators
+ if dd_in.domain_tag is not sym.FRESTR_INTERIOR_FACES:
+ raise ValueError("dd_in must be an interior faces domain")
+ if dd_out != dd_in:
+ raise ValueError("dd_out and dd_in must be identical")
-class FluxOperatorBase(Operator):
- def __init__(self, flux, is_lift=False):
- Operator.__init__(self)
- self.flux = flux
- self.is_lift = is_lift
+ super(OppositeInteriorFaceSwap, self).__init__(dd_in, dd_out)
- def get_flux_or_lift_text(self):
- if self.is_lift:
- return "Lift"
- else:
- return "Flux"
+ mapper_method = intern("map_opposite_interior_face_swap")
- def repr_op(self):
- """Return an equivalent operator with the flux expression set to 0."""
- return type(self)(0, *self.__getinitargs__()[1:])
- def __call__(self, arg):
- # override to suppress apply-operator-to-each-operand
- # behavior from superclass
+class FaceMassOperator(ElementwiseLinearOperator):
+ def __init__(self, dd_in=None, dd_out=None):
+ sym = _sym()
- from grudge.symbolic.primitives import OperatorBinding
- return OperatorBinding(self, arg)
+ if dd_in is None:
+ dd_in = sym.DOFDesc(sym.FRESTR_ALL_FACES, None)
- def __mul__(self, arg):
- from warnings import warn
- warn("Multiplying by a flux operator is deprecated. "
- "Use the less ambiguous parenthesized syntax instead.",
- DeprecationWarning, stacklevel=2)
- return self.__call__(arg)
+ if dd_out is None:
+ dd_out = sym.DOFDesc("vol", None)
+ if not dd_out.is_volume():
+ raise ValueError("dd_out must be a volume domain")
+ if dd_in.domain_tag is not sym.FRESTR_ALL_FACES:
+ raise ValueError("dd_in must be an interior faces domain")
-class QuadratureFluxOperatorBase(FluxOperatorBase):
- pass
-
+ super(FaceMassOperator, self).__init__(dd_in, dd_out)
-class BoundaryFluxOperatorBase(FluxOperatorBase):
- pass
-
-
-class FluxOperator(FluxOperatorBase):
- def __getinitargs__(self):
- return (self.flux, self.is_lift)
-
- mapper_method = intern("map_flux")
-
-
-class BoundaryFluxOperator(BoundaryFluxOperatorBase):
- """
- .. note::
-
- This operator is purely for internal use. It is inserted
- by :class:`grudge.symbolic.mappers.OperatorSpecializer`
- when a :class:`FluxOperator` is applied to a boundary field.
- """
- def __init__(self, flux, boundary_tag, is_lift=False):
- FluxOperatorBase.__init__(self, flux, is_lift)
- self.boundary_tag = boundary_tag
-
- def __getinitargs__(self):
- return (self.flux, self.boundary_tag, self.is_lift)
-
- mapper_method = intern("map_bdry_flux")
-
-
-class QuadratureFluxOperator(QuadratureFluxOperatorBase):
- """
- .. note::
-
- This operator is purely for internal use. It is inserted
- by :class:`grudge.symbolic.mappers.OperatorSpecializer`
- when a :class:`FluxOperator` is applied to a quadrature field.
- """
-
- def __init__(self, flux, quadrature_tag):
- FluxOperatorBase.__init__(self, flux, is_lift=False)
-
- self.quadrature_tag = quadrature_tag
-
- def __getinitargs__(self):
- return (self.flux, self.quadrature_tag)
-
- mapper_method = intern("map_quad_flux")
-
-
-class QuadratureBoundaryFluxOperator(
- QuadratureFluxOperatorBase, BoundaryFluxOperatorBase):
- """
- .. note::
-
- This operator is purely for internal use. It is inserted
- by :class:`grudge.symbolic.mappers.OperatorSpecializer`
- when a :class:`FluxOperator` is applied to a quadrature
- boundary field.
- """
- def __init__(self, flux, quadrature_tag, boundary_tag):
- FluxOperatorBase.__init__(self, flux, is_lift=False)
- self.quadrature_tag = quadrature_tag
- self.boundary_tag = boundary_tag
-
- def __getinitargs__(self):
- return (self.flux, self.quadrature_tag, self.boundary_tag)
-
- mapper_method = intern("map_quad_bdry_flux")
-
-
-class VectorFluxOperator(object):
- """Note that this isn't an actual operator. It's just a placeholder that pops
- out a vector of FluxOperators when applied to an operand.
- """
- def __init__(self, fluxes):
- self.fluxes = fluxes
-
- def __call__(self, arg):
- if isinstance(arg, int) and arg == 0:
- return 0
- from pytools.obj_array import make_obj_array
- from grudge.symbolic.primitives import OperatorBinding
-
- return make_obj_array(
- [OperatorBinding(FluxOperator(f), arg)
- for f in self.fluxes])
-
- def __mul__(self, arg):
- from warnings import warn
- warn("Multiplying by a vector flux operator is deprecated. "
- "Use the less ambiguous parenthesized syntax instead.",
- DeprecationWarning, stacklevel=2)
- return self.__call__(arg)
-
-
-class WholeDomainFluxOperator(pymbolic.primitives.AlgebraicLeaf):
- """Used by the CUDA backend to represent a flux computation on the
- whole domain--interior and boundary.
-
- Unlike other operators, :class:`WholeDomainFluxOperator` instances
- are not bound.
- """
-
- class FluxInfo(Record):
- __slots__ = []
-
- def __repr__(self):
- # override because we want flux_expr in infix
- return "%s(%s)" % (
- self.__class__.__name__,
- ", ".join("%s=%s" % (fld, getattr(self, fld))
- for fld in self.__class__.fields
- if hasattr(self, fld)))
-
- class InteriorInfo(FluxInfo):
- # attributes: flux_expr, field_expr,
-
- @property
- @memoize_method
- def dependencies(self):
- from grudge.symbolic.tools import get_flux_dependencies
- return set(get_flux_dependencies(
- self.flux_expr, self.field_expr))
-
- class BoundaryInfo(FluxInfo):
- # attributes: flux_expr, bpair
-
- @property
- @memoize_method
- def int_dependencies(self):
- from grudge.symbolic.tools import get_flux_dependencies
- return set(get_flux_dependencies(
- self.flux_expr, self.bpair, bdry="int"))
-
- @property
- @memoize_method
- def ext_dependencies(self):
- from grudge.symbolic.tools import get_flux_dependencies
- return set(get_flux_dependencies(
- self.flux_expr, self.bpair, bdry="ext"))
-
- def __init__(self, is_lift, interiors, boundaries,
- quadrature_tag):
- from grudge.symbolic.tools import get_flux_dependencies
-
- self.is_lift = is_lift
-
- self.interiors = tuple(interiors)
- self.boundaries = tuple(boundaries)
- self.quadrature_tag = quadrature_tag
-
- from pytools import set_sum
- interior_deps = set_sum(iflux.dependencies
- for iflux in interiors)
- boundary_int_deps = set_sum(bflux.int_dependencies
- for bflux in boundaries)
- boundary_ext_deps = set_sum(bflux.ext_dependencies
- for bflux in boundaries)
-
- self.interior_deps = list(interior_deps)
- self.boundary_int_deps = list(boundary_int_deps)
- self.boundary_ext_deps = list(boundary_ext_deps)
- self.boundary_deps = list(boundary_int_deps | boundary_ext_deps)
-
- self.dep_to_tag = {}
- for bflux in boundaries:
- for dep in get_flux_dependencies(
- bflux.flux_expr, bflux.bpair, bdry="ext"):
- self.dep_to_tag[dep] = bflux.bpair.tag
-
- def stringifier(self):
- from grudge.symbolic.mappers import StringifyMapper
- return StringifyMapper
-
- def repr_op(self):
- return type(self)(False, [], [], self.quadrature_tag)
-
- @memoize_method
- def rebuild_optemplate(self):
- def generate_summands():
- for i in self.interiors:
- if self.quadrature_tag is None:
- yield FluxOperator(
- i.flux_expr, self.is_lift)(i.field_expr)
- else:
- yield QuadratureFluxOperator(
- i.flux_expr, self.quadrature_tag)(i.field_expr)
- for b in self.boundaries:
- if self.quadrature_tag is None:
- yield BoundaryFluxOperator(
- b.flux_expr, b.bpair.tag, self.is_lift)(b.bpair)
- else:
- yield QuadratureBoundaryFluxOperator(
- b.flux_expr, self.quadrature_tag,
- b.bpair.tag)(b.bpair)
-
- from pymbolic.primitives import flattened_sum
- return flattened_sum(generate_summands())
-
- # infrastructure interaction
- def get_hash(self):
- return hash((self.__class__, self.rebuild_optemplate()))
-
- def is_equal(self, other):
- return (other.__class__ == WholeDomainFluxOperator
- and self.rebuild_optemplate() == other.rebuild_optemplate())
-
- def __getinitargs__(self):
- return (self.is_lift, self.interiors, self.boundaries,
- self.quadrature_tag)
-
- mapper_method = intern("map_whole_domain_flux")
+ mapper_method = intern("map_face_mass_operator")
# }}}
# {{{ convenience functions for operator creation
-def get_flux_operator(flux):
- """Return a flux operator that can be multiplied with
- a volume field to obtain the interior fluxes
- or with a :class:`BoundaryPair` to obtain the lifted boundary
- flux.
- """
- from pytools.obj_array import is_obj_array
- from grudge.symbolic.operators import VectorFluxOperator, FluxOperator
-
- if is_obj_array(flux):
- return VectorFluxOperator(flux)
- else:
- return FluxOperator(flux)
-
-
def nabla(dim):
from pytools.obj_array import make_obj_array
return make_obj_array(
- [DifferentiationOperator(i) for i in range(dim)])
+ [DiffOperator(i) for i in range(dim)])
def minv_stiffness_t(dim):
@@ -760,19 +431,25 @@ def stiffness_t(dim):
[StiffnessTOperator(i) for i in range(dim)])
-def integral(arg):
+def integral(arg, dd=None):
+ if dd is None:
+ dd = _sym().DD_VOLUME
+
from grudge import sym
- return sym.NodalSum()(sym.MassOperator()(sym.Ones())*arg)
+ return sym.NodalSum(dd)(sym.MassOperator()(sym.Ones(dd))*arg)
-def norm(p, arg):
+def norm(p, arg, dd=None):
"""
:arg arg: is assumed to be a vector, i.e. have shape ``(n,)``.
"""
import grudge.symbolic as sym
+ if dd is None:
+ dd = _sym().DD_VOLUME
+
if p == 2:
- comp_norm_squared = sym.NodalSum()(
+ comp_norm_squared = sym.NodalSum(dd_in=dd)(
sym.CFunction("fabs")(
arg * sym.MassOperator()(arg)))
return sym.CFunction("sqrt")(sum(comp_norm_squared))
@@ -783,7 +460,7 @@ def norm(p, arg):
return reduce(Max, comp_norm)
else:
- return sum(sym.NodalSum()(sym.CFunction("fabs")(arg)**p))**(1/p)
+ return sum(sym.NodalSum(dd_in=dd)(sym.CFunction("fabs")(arg)**p))**(1/p)
# }}}
diff --git a/grudge/symbolic/primitives.py b/grudge/symbolic/primitives.py
index 5be9bfda..0411e888 100644
--- a/grudge/symbolic/primitives.py
+++ b/grudge/symbolic/primitives.py
@@ -29,6 +29,8 @@ from six.moves import range, intern
import numpy as np
import pymbolic.primitives
from meshmode.mesh import BTAG_ALL, BTAG_NONE # noqa
+from meshmode.discretization.connection import ( # noqa
+ FRESTR_ALL_FACES, FRESTR_INTERIOR_FACES)
from pymbolic.primitives import ( # noqa
cse_scope as cse_scope_base,
@@ -37,24 +39,40 @@ from pymbolic.geometric_algebra import MultiVector
from pytools.obj_array import join_fields, make_obj_array # noqa
-class LeafBase(pymbolic.primitives.AlgebraicLeaf):
+class ExpressionBase(pymbolic.primitives.Expression):
def stringifier(self):
from grudge.symbolic.mappers import StringifyMapper
return StringifyMapper
-class VTAG_ALL:
- """This is used in a 'where' field to denote the volume discretization."""
- pass
+def _sym():
+ # A hack to make referring to grudge.sym possible without
+ # circular imports and tons of typing.
+
+ from grudge import sym
+ return sym
__doc__ = """
+
+.. currentmodule:: grudge.sym
+
+DOF description
+^^^^^^^^^^^^^^^
+
+.. autoclass:: DTAG_SCALAR
+.. autoclass:: DTAG_VOLUME_ALL
+.. autoclass:: QTAG_NONE
+.. autoclass:: DOFDesc
+.. data:: DD_SCALAR
+.. data:: DD_VOLUME
+
Symbols
^^^^^^^
-.. autoclass:: Field
-.. autoclass:: make_sym_vector
+.. autoclass:: Variable
.. autoclass:: make_sym_array
+.. autoclass:: make_sym_mv
.. autoclass:: ScalarParameter
.. autoclass:: CFunction
@@ -82,33 +100,205 @@ Geometry data
"""
+# {{{ DOF description
+
+class DTAG_SCALAR:
+ pass
+
+
+class DTAG_VOLUME_ALL:
+ pass
+
+
+class QTAG_NONE:
+ pass
+
+
+class DOFDesc(object):
+ """Describes the meaning of degrees of freedom.
+
+ .. attribute:: domain_tag
+ .. attribute:: quadrature_tag
+ .. automethod:: is_scalar
+ .. automethod:: is_discretized
+ .. automethod:: is_volume
+ .. automethod:: is_boundary
+ .. automethod:: uses_quadrature
+ .. automethod:: with_qtag
+ .. automethod:: with_dtag
+ .. automethod:: __eq__
+ .. automethod:: __ne__
+ .. automethod:: __hash__
+ """
+
+ def __init__(self, domain_tag, quadrature_tag=None):
+ """
+ :arg domain_tag: One of the following:
+ :class:`grudge.sym.DTAG_SCALAR` (or the string ``"scalar"``),
+ :class:`grudge.sym.DTAG_VOLUME_ALL` (or the string ``"vol"``)
+ for the default volume discretization,
+ :class:`meshmode.discretization.connection.FRESTR_ALL_FACES`
+ (or the string ``"all_faces"``),
+ or
+ :class:`meshmode.discretization.connection.FRESTR_INTERIOR_FACES`
+ (or the string ``"int_faces"``),
+ or *None* to indicate that the geometry is not yet known.
+
+ :arg quadrature_tag:
+ *None* to indicate that the quadrature grid is not known,or
+ :class:`QTAG_NONE` to indicate the use of the basic discretization
+ grid, or a string to indicate the use of the thus-tagged quadratue
+ grid.
+ """
+ if domain_tag == "scalar":
+ domain_tag = DTAG_SCALAR
+ elif domain_tag is DTAG_SCALAR:
+ domain_tag = DTAG_SCALAR
+ elif domain_tag == "vol":
+ domain_tag = DTAG_VOLUME_ALL
+ elif domain_tag is DTAG_VOLUME_ALL:
+ pass
+ elif domain_tag == "all_faces":
+ domain_tag = FRESTR_ALL_FACES
+ elif domain_tag is FRESTR_ALL_FACES:
+ pass
+ elif domain_tag == "int_faces":
+ domain_tag = FRESTR_INTERIOR_FACES
+ elif domain_tag is FRESTR_INTERIOR_FACES:
+ pass
+ elif domain_tag is None:
+ pass
+
+ if domain_tag is DTAG_SCALAR and quadrature_tag is not None:
+ raise ValueError("cannot have nontrivial quadrature tag on scalar")
+
+ self.domain_tag = domain_tag
+ self.quadrature_tag = quadrature_tag
+
+ def is_scalar(self):
+ return self.domain_tag is DTAG_SCALAR
+
+ def is_discretized(self):
+ return not self.is_scalar()
+
+ def is_volume(self):
+ return self.domain_tag is DTAG_VOLUME_ALL
+
+ def is_boundary(self):
+ return not self.is_volume()
+
+ def uses_quadrature(self):
+ if self.quadrature_tag is None:
+ return False
+ if self.quadrature_tag is QTAG_NONE:
+ return False
+
+ return True
+
+ def with_qtag(self, qtag):
+ return type(self)(domain_tag=self.domain_tag, quadrature_tag=qtag)
+
+ def with_dtag(self, dtag):
+ return type(self)(domain_tag=dtag, quadrature_tag=self.quadrature_tag)
+
+ def __eq__(self, other):
+ return (type(self) == type(other)
+ and self.domain_tag == other.domain_tag
+ and self.quadrature_tag == other.quadrature_tag)
+
+ def __ne__(self, other):
+ return not self.__eq__(other)
+
+ def __hash__(self):
+ return hash((type(self), self.domain_tag, self.quadrature_tag))
+
+ def __repr__(self):
+ def fmt(s):
+ if isinstance(s, type):
+ return s.__name__
+ else:
+ return repr(s)
+
+ return "DOFDesc(%s, %s)" % (fmt(self.domain_tag), fmt(self.quadrature_tag))
+
+
+DD_SCALAR = DOFDesc(DTAG_SCALAR, None)
+
+DD_VOLUME = DOFDesc(DTAG_VOLUME_ALL, None)
+
+
+def as_dofdesc(dd):
+ if isinstance(dd, DOFDesc):
+ return dd
+ else:
+ return DOFDesc(dd, None)
+
+# }}}
+
+
+# {{{ has-dof-desc mix-in
+
+class HasDOFDesc(object):
+ """
+ .. attribute:: dd
+
+ an instance of :class:`grudge.sym.DOFDesc` describing the
+ discretization on which this property is given.
+ """
+
+ def __init__(self, dd):
+ self.dd = dd
+
+ def __getinitargs__(self):
+ return (self.dd,)
+
+ def with_dd(self, dd):
+ """Return a copy of *self*, modified to the given DOF descriptor.
+ """
+ return type(self)(
+ *self.__getinitargs__()[:-1],
+ dd=dd or self.dd)
+
+# }}}
+
+
# {{{ variables
class cse_scope(cse_scope_base): # noqa
DISCRETIZATION = "grudge_discretization"
-Field = pymbolic.primitives.Variable
+class Variable(HasDOFDesc, ExpressionBase, pymbolic.primitives.Variable):
+ def __init__(self, name, dd=None):
+ if dd is None:
+ dd = DD_VOLUME
-make_sym_vector = pymbolic.primitives.make_sym_vector
-make_sym_array = pymbolic.primitives.make_sym_array
+ HasDOFDesc.__init__(self, dd)
+ pymbolic.primitives.Variable.__init__(self, name)
+ def __getinitargs__(self):
+ return (self.name, self.dd,)
-def make_field(var_or_string):
- if not isinstance(var_or_string, pymbolic.primitives.Expression):
- return Field(var_or_string)
- else:
- return var_or_string
+ mapper_method = "map_grudge_variable"
+var = Variable
-class ScalarParameter(pymbolic.primitives.Variable):
- """A placeholder for a user-supplied scalar variable."""
- def stringifier(self):
- from grudge.symbolic.mappers import StringifyMapper
- return StringifyMapper
+class ScalarVariable(Variable):
+ def __init__(self, name):
+ super(ScalarVariable, self).__init__(name, DD_SCALAR)
+
- mapper_method = intern("map_scalar_parameter")
+def make_sym_array(name, shape, dd=None):
+ def var_factory(name):
+ return Variable(name, dd)
+
+ return pymbolic.primitives.make_sym_array(name, shape, var_factory)
+
+
+def make_sym_mv(name, dim, var_factory=None):
+ return MultiVector(
+ make_sym_array(name, dim, var_factory))
class CFunction(pymbolic.primitives.Variable):
@@ -139,7 +329,7 @@ cos = CFunction("cos")
# {{{ technical helpers
-class OperatorBinding(LeafBase):
+class OperatorBinding(ExpressionBase):
def __init__(self, op, field):
self.op = op
self.field = field
@@ -179,99 +369,47 @@ class PrioritizedSubexpression(pymbolic.primitives.CommonSubexpression):
def get_extra_properties(self):
return {"priority": self.priority}
-
-class BoundaryPair(LeafBase):
- """Represents a pairing of a volume and a boundary field, used for the
- application of boundary fluxes.
- """
-
- def __init__(self, field, bfield, tag=BTAG_ALL):
- self.field = field
- self.bfield = bfield
- self.tag = tag
-
- mapper_method = intern("map_boundary_pair")
-
- def __getinitargs__(self):
- return (self.field, self.bfield, self.tag)
-
- def get_hash(self):
- from pytools.obj_array import obj_array_to_hashable
-
- return hash((self.__class__,
- obj_array_to_hashable(self.field),
- obj_array_to_hashable(self.bfield),
- self.tag))
-
- def is_equal(self, other):
- from pytools.obj_array import obj_array_equal
- return (self.__class__ == other.__class__
- and obj_array_equal(other.field, self.field)
- and obj_array_equal(other.bfield, self.bfield)
- and other.tag == self.tag)
-
# }}}
-class Ones(LeafBase):
- def __init__(self, quadrature_tag=None, where=None):
- self.where = where
- self.quadrature_tag = quadrature_tag
-
- def __getinitargs__(self):
- return (self.where, self.quadrature_tag,)
-
+class Ones(ExpressionBase, HasDOFDesc):
mapper_method = intern("map_ones")
# {{{ geometry data
-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.
-
- .. attribute:: quadrature_tag
-
- quadrature tag for the grid on
- which this geometric factor is needed, or None for
- nodal representation.
- """
-
- def __init__(self, quadrature_tag, where=None):
- self.quadrature_tag = quadrature_tag
- self.where = where
-
- def __getinitargs__(self):
- return (self.quadrature_tag, self.where)
+class DiscretizationProperty(ExpressionBase, HasDOFDesc):
+ pass
class NodeCoordinateComponent(DiscretizationProperty):
+ def __init__(self, axis, dd=None):
+ if dd is None:
+ dd = DD_VOLUME
+
+ if not dd.is_discretized():
+ raise ValueError("dd must be a discretization for "
+ "NodeCoordinateComponent")
- def __init__(self, axis, quadrature_tag=None, where=None):
- super(NodeCoordinateComponent, self).__init__(quadrature_tag, where)
+ super(NodeCoordinateComponent, self).__init__(dd)
self.axis = axis
def __getinitargs__(self):
- return (self.axis, self.quadrature_tag)
+ return (self.axis, self.dd)
mapper_method = intern("map_node_coordinate_component")
-def nodes(ambient_dim, quadrature_tag=None, where=None):
- return np.array([NodeCoordinateComponent(i, quadrature_tag, where)
+def nodes(ambient_dim, dd=None):
+ return np.array([NodeCoordinateComponent(i, dd)
for i in range(ambient_dim)], dtype=object)
-def mv_nodes(ambient_dim, quadrature_tag=None, where=None):
- return MultiVector(
- nodes(ambient_dim, quadrature_tag=quadrature_tag, where=where))
+def mv_nodes(ambient_dim, dd=None):
+ return MultiVector(nodes(ambient_dim, dd))
-def forward_metric_derivative(xyz_axis, rst_axis, where=None,
- quadrature_tag=None):
+def forward_metric_derivative(xyz_axis, rst_axis, dd=None):
r"""
Pointwise metric derivatives representing
@@ -279,15 +417,17 @@ def forward_metric_derivative(xyz_axis, rst_axis, where=None,
\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)
+ if dd is None:
+ dd = DD_VOLUME
+
+ inner_dd = dd.with_qtag(QTAG_NONE)
- result = diff_op(NodeCoordinateComponent(xyz_axis, where=where))
+ diff_op = _sym().RefDiffOperator(rst_axis, inner_dd)
- if quadrature_tag is not None:
- result = QuadratureGridUpsampler(quadrature_tag, where)(result)
+ result = diff_op(NodeCoordinateComponent(xyz_axis, inner_dd))
+
+ if dd.uses_quadrature():
+ result = _sym().interp(inner_dd, dd)(result)
return cse(
result,
@@ -295,35 +435,29 @@ def forward_metric_derivative(xyz_axis, rst_axis, where=None,
scope=cse_scope.DISCRETIZATION)
-def forward_metric_derivative_vector(ambient_dim, rst_axis, where=None,
- quadrature_tag=None):
+def forward_metric_derivative_vector(ambient_dim, rst_axis, dd=None):
return make_obj_array([
- forward_metric_derivative(
- i, rst_axis, where=where, quadrature_tag=quadrature_tag)
+ forward_metric_derivative(i, rst_axis, dd=dd)
for i in range(ambient_dim)])
-def forward_metric_derivative_mv(ambient_dim, rst_axis, where=None,
- quadrature_tag=None):
+def forward_metric_derivative_mv(ambient_dim, rst_axis, dd=None):
return MultiVector(
- forward_metric_derivative_vector(
- ambient_dim, rst_axis, where=where, quadrature_tag=quadrature_tag))
+ forward_metric_derivative_vector(ambient_dim, rst_axis, dd=dd))
-def parametrization_derivative(ambient_dim, dim=None, where=None,
- quadrature_tag=None):
+def parametrization_derivative(ambient_dim, dim=None, dd=None):
if dim is None:
dim = ambient_dim
from pytools import product
return product(
- forward_metric_derivative_mv(
- ambient_dim, rst_axis, where, quadrature_tag)
+ forward_metric_derivative_mv(ambient_dim, rst_axis, dd)
for rst_axis in range(dim))
def inverse_metric_derivative(rst_axis, xyz_axis, ambient_dim, dim=None,
- where=None, quadrature_tag=None):
+ dd=None):
if dim is None:
dim = ambient_dim
@@ -332,8 +466,7 @@ def inverse_metric_derivative(rst_axis, xyz_axis, ambient_dim, dim=None,
"the derivative matrix is not square")
par_vecs = [
- forward_metric_derivative_mv(
- ambient_dim, rst, where, quadrature_tag)
+ forward_metric_derivative_mv(ambient_dim, rst, dd)
for rst in range(dim)]
# Yay Cramer's rule! (o rly?)
@@ -352,7 +485,7 @@ def inverse_metric_derivative(rst_axis, xyz_axis, ambient_dim, dim=None,
volume_pseudoscalar_inv = cse(outerprod(
forward_metric_derivative_mv(
- ambient_dim, rst_axis, where, quadrature_tag)
+ ambient_dim, rst_axis, dd=dd)
for rst_axis in range(dim)).inv())
return cse(
@@ -362,20 +495,17 @@ def inverse_metric_derivative(rst_axis, xyz_axis, ambient_dim, dim=None,
scope=cse_scope.DISCRETIZATION)
-def forward_metric_derivative_mat(ambient_dim, dim=None,
- where=None, quadrature_tag=None):
+def forward_metric_derivative_mat(ambient_dim, dim=None, dd=None):
if dim is None:
dim = ambient_dim
result = np.zeros((ambient_dim, dim), dtype=np.object)
for j in range(dim):
- result[:, j] = forward_metric_derivative_vector(ambient_dim, j,
- where=where, quadrature_tag=quadrature_tag)
+ result[:, j] = forward_metric_derivative_vector(ambient_dim, j, dd=dd)
return result
-def inverse_metric_derivative_mat(ambient_dim, dim=None,
- where=None, quadrature_tag=None):
+def inverse_metric_derivative_mat(ambient_dim, dim=None, dd=None):
if dim is None:
dim = ambient_dim
@@ -383,34 +513,37 @@ def inverse_metric_derivative_mat(ambient_dim, dim=None,
for i in range(dim):
for j in range(ambient_dim):
result[i, j] = inverse_metric_derivative(
- i, j, ambient_dim, dim,
- where=where, quadrature_tag=quadrature_tag)
+ i, j, ambient_dim, dim, dd=dd)
return result
-def pseudoscalar(ambient_dim, dim=None, where=None, quadrature_tag=None):
+def pseudoscalar(ambient_dim, dim=None, dd=None):
if dim is None:
dim = ambient_dim
return cse(
- parametrization_derivative(ambient_dim, dim, where=where,
- quadrature_tag=quadrature_tag)
+ parametrization_derivative(ambient_dim, dim, dd=dd)
.project_max_grade(),
"pseudoscalar", cse_scope.DISCRETIZATION)
-def area_element(ambient_dim, dim=None, where=None, quadrature_tag=None):
+def area_element(ambient_dim, dim=None, dd=None):
return cse(
sqrt(
- pseudoscalar(ambient_dim, dim, where, quadrature_tag=quadrature_tag)
+ pseudoscalar(ambient_dim, dim, dd=dd)
.norm_squared()),
"area_element", cse_scope.DISCRETIZATION)
-def mv_normal(tag, ambient_dim, dim=None, quadrature_tag=None):
+def mv_normal(dd, ambient_dim, dim=None):
"""Exterior unit normal as a MultiVector."""
+ dd = as_dofdesc(dd)
+
+ if not dd.is_boundary():
+ raise ValueError("may only request normals on boundaries")
+
if dim is None:
dim = ambient_dim - 1
@@ -418,16 +551,95 @@ def mv_normal(tag, ambient_dim, dim=None, quadrature_tag=None):
# exterior normals for positively oriented curves.
pder = (
- pseudoscalar(ambient_dim, dim, tag, quadrature_tag=quadrature_tag)
+ pseudoscalar(ambient_dim, dim, dd=dd)
/
- area_element(ambient_dim, dim, tag, quadrature_tag=quadrature_tag))
+ area_element(ambient_dim, dim, dd=dd))
return cse(pder.I | pder, "normal",
cse_scope.DISCRETIZATION)
-def normal(tag, ambient_dim, dim=None, quadrature_tag=None):
- return mv_normal(tag, ambient_dim, dim,
- quadrature_tag=quadrature_tag).as_vector()
+def normal(dd, ambient_dim, dim=None, quadrature_tag=None):
+ return mv_normal(dd, ambient_dim, dim).as_vector()
+
+# }}}
+
+
+# {{{ flux parir
+
+class FluxPair:
+ """
+ .. attribute:: dd
+
+ an instance of :class:`grudge.sym.DOFDesc` describing the
+ discretization on which :attr:`interior` and :attr:`exterior`
+ live.
+
+ .. attribute:: interior
+
+ an expression representing the interior value to
+ be used for the flux.
+
+ .. attribute:: exterior
+
+ an expression representing the exterior value to
+ be used for the flux.
+ """
+ def __init__(self, dd, interior, exterior):
+ """
+ """
+
+ self.dd = as_dofdesc(dd)
+ self.interior = interior
+ self.exterior = exterior
+
+ def __getitem__(self, index):
+ return FluxPair(
+ self.dd,
+ self.exterior[index],
+ self.interior[index])
+
+ @property
+ def int(self):
+ return self.interior
+
+ @property
+ def ext(self):
+ return self.exterior
+
+ @property
+ def avg(self):
+ return 0.5*(self.int + self.ext)
+
+
+def int_fpair(expression):
+ i = cse(_sym().interp("vol", "int_faces")(expression))
+ e = cse(_sym().OppositeInteriorFaceSwap()(i))
+ return FluxPair("int_faces", i, e)
+
+
+def bdry_fpair(dd, interior, exterior):
+ """
+ :arg interior: an expression that already lives on the boundary
+ representing the interior value to be used
+ for the flux.
+ :arg exterior: an expression that already lives on the boundary
+ representing the exterior value to be used
+ for the flux.
+ """
+ return FluxPair(dd, interior, exterior)
+
+
+def bv_fpair(dd, interior, exterior):
+ """
+ :arg interior: an expression that lives in the volume
+ and will be restricted to the boundary identified by
+ *tag* before use.
+ :arg exterior: an expression that already lives on the boundary
+ representing the exterior value to be used
+ for the flux.
+ """
+ interior = _sym().cse(_sym().interp("vol", dd)(interior))
+ return FluxPair(dd, interior, exterior)
# }}}
diff --git a/grudge/symbolic/tools.py b/grudge/symbolic/tools.py
index 524b7cef..55a654d6 100644
--- a/grudge/symbolic/tools.py
+++ b/grudge/symbolic/tools.py
@@ -33,48 +33,10 @@ import numpy as np
def is_scalar(expr):
from grudge import sym
- return isinstance(expr, (int, float, complex, sym.ScalarParameter))
+ if isinstance(expr, sym.Variable) and expr.dd.domain_tag is sym.DTAG_SCALAR:
+ return True
-
-def get_flux_dependencies(flux, field, bdry="all"):
- from grudge.symbolic.flux.mappers import FluxDependencyMapper
- from grudge.symbolic.flux.primitives import FieldComponent
- in_fields = list(FluxDependencyMapper(
- include_calls=False)(flux))
-
- # check that all in_fields are FieldComponent instances
- assert not [in_field
- for in_field in in_fields
- if not isinstance(in_field, FieldComponent)]
-
- def maybe_index(fld, index):
- from pytools.obj_array import is_obj_array
- if is_obj_array(fld):
- return fld[inf.index]
- else:
- return fld
-
- from grudge.tools import is_zero
- from grudge import sym
- if isinstance(field, sym.BoundaryPair):
- for inf in in_fields:
- if inf.is_interior:
- if bdry in ["all", "int"]:
- value = maybe_index(field.field, inf.index)
-
- if not is_zero(value):
- yield value
- else:
- if bdry in ["all", "ext"]:
- value = maybe_index(field.bfield, inf.index)
-
- if not is_zero(value):
- yield value
- else:
- for inf in in_fields:
- value = maybe_index(field, inf.index)
- if not is_zero(value):
- yield value
+ return isinstance(expr, (int, float, complex))
def split_sym_operator_for_multirate(state_vector, sym_operator,
@@ -169,22 +131,10 @@ def pretty(sym_operator):
splitter = "="*75 + "\n"
- bc_strs = stringify_mapper.get_bc_strings()
- if bc_strs:
- result = "\n".join(bc_strs)+"\n"+splitter+result
-
cse_strs = stringify_mapper.get_cse_strings()
if cse_strs:
result = "\n".join(cse_strs)+"\n"+splitter+result
- flux_strs = stringify_mapper.get_flux_strings()
- if flux_strs:
- result = "\n".join(flux_strs)+"\n"+splitter+result
-
- flux_cses = stringify_mapper.flux_stringify_mapper.get_cse_strings()
- if flux_cses:
- result = "\n".join("flux "+fs for fs in flux_cses)+"\n\n"+result
-
return result
# }}}
diff --git a/test/test_grudge.py b/test/test_grudge.py
index 3754f49f..c137fe26 100644
--- a/test/test_grudge.py
+++ b/test/test_grudge.py
@@ -102,13 +102,13 @@ def test_1d_mass_mat_trig(ctx_getter):
x = sym.nodes(1)
f = bind(discr, sym.cos(x[0])**2)(queue)
- ones = bind(discr, sym.Ones())(queue)
+ ones = bind(discr, sym.Ones(sym.DD_VOLUME))(queue)
- mass_op = bind(discr, sym.MassOperator()(sym.Field("f")))
+ mass_op = bind(discr, sym.MassOperator()(sym.var("f")))
num_integral_1 = np.dot(ones.get(), mass_op(queue, f=f).get())
num_integral_2 = np.dot(f.get(), mass_op(queue, f=ones).get())
- num_integral_3 = bind(discr, sym.integral(sym.Field("f")))(queue, f=f).get()
+ num_integral_3 = bind(discr, sym.integral(sym.var("f")))(queue, f=f).get()
true_integral = 13*np.pi/2
err_1 = abs(num_integral_1-true_integral)
@@ -148,7 +148,7 @@ def test_tri_diff_mat(ctx_getter, dim, order=4):
f = bind(discr, sym.sin(3*x[axis]))(queue)
df = bind(discr, 3*sym.cos(3*x[axis]))(queue)
- sym_op = nabla[axis](sym.Field("f"))
+ sym_op = nabla[axis](sym.var("f"))
bound_op = bind(discr, sym_op)
df_num = bound_op(queue, f=f)
--
GitLab