diff --git a/grudge/discretization.py b/grudge/discretization.py
index 300b5fbc362b4d84ff940a72efd0d686d99dc135..6a7a4fb453e8ff0ae60295eb3602ee80416c266b 100644
--- a/grudge/discretization.py
+++ b/grudge/discretization.py
@@ -83,11 +83,13 @@ class Discretization(object):
def mesh(self):
return self.volume_discr.mesh
- def empty(self, queue=None, dtype=None, extra_dims=None):
- return self.volume_discr.empty(queue, dtype, extra_dims=None)
+ def empty(self, queue=None, dtype=None, extra_dims=None, allocator=None):
+ return self.volume_discr.empty(queue, dtype, extra_dims=None,
+ allocator=allocator)
- def zeros(self, queue, dtype=None, extra_dims=None):
- return self.volume_discr.zeros(queue, dtype, extra_dims=None)
+ def zeros(self, queue, dtype=None, extra_dims=None, allocator=None):
+ return self.volume_discr.zeros(queue, dtype, extra_dims=None,
+ allocator=allocator)
def is_volume_where(self, where):
from grudge import sym
diff --git a/grudge/execution.py b/grudge/execution.py
index 4f0293751d79eb91d0dc8837f71140c92b5841f2..2cf24b316cbfdf80b1859a2f00ff9123d213822e 100644
--- a/grudge/execution.py
+++ b/grudge/execution.py
@@ -24,6 +24,11 @@ THE SOFTWARE.
import six
import numpy as np
+import loopy as lp
+import pyopencl as cl
+import pyopencl.array # noqa
+from pytools import memoize_in
+
import grudge.symbolic.mappers as mappers
from grudge import sym
@@ -45,17 +50,17 @@ class ExecutionMapper(mappers.Evaluator,
# {{{ expression mappings -------------------------------------------------
def map_ones(self, expr):
- # FIXME
if expr.quadrature_tag is not None:
+ # FIXME
raise NotImplementedError("ones on quad. grids")
- result = self.discr.empty(self.queue)
+ result = self.discr.empty(self.queue, allocator=self.bound_op.allocator)
result.fill(1)
return result
def map_node_coordinate_component(self, expr):
- # FIXME
if expr.quadrature_tag is not None:
+ # FIXME
raise NotImplementedError("node coordinate components on quad. grids")
if self.discr.is_volume_where(expr.where):
@@ -84,32 +89,25 @@ class ExecutionMapper(mappers.Evaluator,
return func(*[self.rec(p) for p in expr.parameters])
def map_nodal_sum(self, op, field_expr):
- return np.sum(self.rec(field_expr))
+ return cl.array.sum(self.rec(field_expr))
def map_nodal_max(self, op, field_expr):
- return np.max(self.rec(field_expr))
+ return cl.array.max(self.rec(field_expr))
def map_nodal_min(self, op, field_expr):
- return np.min(self.rec(field_expr))
+ return cl.array.min(self.rec(field_expr))
def map_if(self, expr):
bool_crit = self.rec(expr.condition)
+
then = self.rec(expr.then)
else_ = self.rec(expr.else_)
- true_indices = np.nonzero(bool_crit)
- false_indices = np.nonzero(~bool_crit)
+ result = cl.array.empty_like(then, queue=self.queue,
+ allocator=self.bound_op.allocator)
+ cl.array.if_positive(bool_crit, then, else_, out=result,
+ queue=self.queue)
- result = self.discr.volume_empty(
- kind=self.discr.compute_kind)
-
- if isinstance(then, np.ndarray):
- then = then[true_indices]
- if isinstance(else_, np.ndarray):
- else_ = else_[false_indices]
-
- result[true_indices] = then
- result[false_indices] = else_
return result
def map_ref_diff_base(self, op, field_expr):
@@ -123,28 +121,43 @@ class ExecutionMapper(mappers.Evaluator,
if is_zero(field):
return 0
- for eg in self.discr.element_groups:
+ @memoize_in(self.bound_op, "elwise_linear_knl")
+ def knl():
+ knl = lp.make_kernel(
+ """{[k,i,j]:
+ 0<=k<nelements and
+ 0<=i,j<ndiscr_nodes}""",
+ "result[k,i] = sum(j, diff_mat[i, j] * vec[k, j])",
+ default_offset=lp.auto, name="diff")
+
+ knl = lp.split_iname(knl, "i", 16, inner_tag="l.0")
+ return lp.tag_inames(knl, dict(k="g.0"))
+
+ discr = self.discr.volume_discr
+
+ # FIXME: This shouldn't really assume that it's dealing with a volume
+ # input. What about quadrature? What about boundaries?
+ result = discr.empty(
+ queue=self.queue,
+ dtype=field.dtype, allocator=self.bound_op.allocator)
+
+ for grp in discr.groups:
+ cache_key = grp, op, field.dtype
try:
- matrix, coeffs = self.bound_op.elwise_linear_cache[eg, op, field.dtype]
+ matrix = self.bound_op.elwise_linear_cache[cache_key]
except KeyError:
- matrix = np.asarray(op.matrix(eg), dtype=field.dtype)
- coeffs = op.coefficients(eg)
- self.elwise_linear_cache[eg, op, field.dtype] = matrix, coeffs
+ matrix = (
+ cl.array.to_device(
+ self.queue,
+ np.asarray(op.matrix(grp), dtype=field.dtype))
+ .with_queue(None))
- from grudge._internal import (
- perform_elwise_scaled_operator,
- perform_elwise_operator)
+ self.bound_op.elwise_linear_cache[cache_key] = matrix
- if coeffs is None:
- perform_elwise_operator(eg.ranges, eg.ranges,
- matrix, field, out)
- else:
- perform_elwise_scaled_operator(eg.ranges, eg.ranges,
- coeffs, matrix, field, out)
+ knl()(self.queue, diff_mat=matrix, result=grp.view(result),
+ vec=grp.view(field))
- out = self.discr.volume_zeros()
- self.bound_op.do_elementwise_linear(op, field, out)
- return out
+ return result
def map_ref_quad_mass(self, op, field_expr):
field = self.rec(field_expr)
@@ -401,6 +414,10 @@ class ExecutionMapper(mappers.Evaluator,
# FIXME
raise NotImplementedError()
+ # FIXME: There's no real reason why differentiation is special,
+ # 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:
@@ -422,11 +439,12 @@ class ExecutionMapper(mappers.Evaluator,
# {{{ bound operator
class BoundOperator(object):
- def __init__(self, discr, code, debug_flags):
+ def __init__(self, discr, code, debug_flags, allocator=None):
self.discr = discr
self.code = code
self.elwise_linear_cache = {}
self.debug_flags = debug_flags
+ self.allocator = allocator
def __call__(self, queue, **context):
import pyopencl.array as cl_array
@@ -527,7 +545,7 @@ def process_sym_operator(sym_operator, post_bind_mapper=None,
def bind(discr, sym_operator, post_bind_mapper=lambda x: x, type_hints={},
- debug_flags=set()):
+ debug_flags=set(), allocator=None):
# from grudge.symbolic.mappers import QuadratureUpsamplerRemover
# sym_operator = QuadratureUpsamplerRemover(self.quad_min_degrees)(
# sym_operator)
@@ -549,9 +567,11 @@ def bind(discr, sym_operator, post_bind_mapper=lambda x: x, type_hints={},
type_hints=type_hints)
from grudge.symbolic.compiler import OperatorCompiler
- code = OperatorCompiler()(sym_operator, type_hints)
+ code = OperatorCompiler()(sym_operator,
+ type_hints=type_hints)
- bound_op = BoundOperator(discr, code, type_hints)
+ bound_op = BoundOperator(discr, code,
+ debug_flags=debug_flags, allocator=allocator)
if "dump_op_code" in debug_flags:
from grudge.tools import open_unique_debug_file
diff --git a/grudge/symbolic/compiler.py b/grudge/symbolic/compiler.py
index 466475229bdc9caf189431c9ace585dd2a6a5de9..55c5fa6d598c40c7bed6a5751e02b7944add9c62 100644
--- a/grudge/symbolic/compiler.py
+++ b/grudge/symbolic/compiler.py
@@ -594,7 +594,7 @@ class Code(object):
self.last_schedule = None
self.static_schedule_attempts -= 1
- from grudge.tools import with_object_array_or_scalar
+ from pytools.obj_array import with_object_array_or_scalar
return with_object_array_or_scalar(exec_mapper, self.result)
# }}}
diff --git a/grudge/symbolic/operators.py b/grudge/symbolic/operators.py
index c63d2f8cdece59e1a17aad20e6f4bbd001ce9efd..c5375c4f68eb0dab3e4f0aad1fa840a835e2ee64 100644
--- a/grudge/symbolic/operators.py
+++ b/grudge/symbolic/operators.py
@@ -226,9 +226,6 @@ class ElementwiseLinearOperator(Operator):
def matrix(self, element_group):
raise NotImplementedError
- def coefficients(self, element_group):
- return None
-
mapper_method = intern("map_elementwise_linear")
@@ -348,31 +345,15 @@ class VandermondeOperator(ElementwiseLinearOperator):
# {{{ mass operators
-class MassOperatorBase(ElementwiseLinearOperator):
+class MassOperatorBase(Operator):
pass
class MassOperator(MassOperatorBase):
- @staticmethod
- def matrix(element_group):
- return element_group.mass_matrix
-
- @staticmethod
- def coefficients(element_group):
- return element_group.jacobians
-
mapper_method = intern("map_mass")
class InverseMassOperator(MassOperatorBase):
- @staticmethod
- def matrix(element_group):
- return element_group.inverse_mass_matrix
-
- @staticmethod
- def coefficients(element_group):
- return element_group.inverse_jacobians
-
mapper_method = intern("map_inverse_mass")
@@ -423,11 +404,10 @@ class ReferenceMassOperatorBase(MassOperatorBase):
class ReferenceMassOperator(ReferenceMassOperatorBase):
@staticmethod
def matrix(element_group):
- return element_group.mass_matrix
-
- @staticmethod
- def coefficients(element_group):
- return None
+ import modepy as mp
+ return mp.mass_matrix(
+ element_group.basis(),
+ element_group.unit_nodes)
mapper_method = intern("map_ref_mass")
@@ -435,11 +415,10 @@ class ReferenceMassOperator(ReferenceMassOperatorBase):
class ReferenceInverseMassOperator(ReferenceMassOperatorBase):
@staticmethod
def matrix(element_group):
- return element_group.inverse_mass_matrix
-
- @staticmethod
- def coefficients(element_group):
- return None
+ import modepy as mp
+ return mp.inverse_mass_matrix(
+ element_group.basis(),
+ element_group.unit_nodes)
mapper_method = intern("map_ref_inverse_mass")
@@ -782,7 +761,7 @@ def stiffness_t(dim):
def integral(arg):
- import grudge.symbolic as sym
+ from grudge import sym
return sym.NodalSum()(sym.MassOperator()(sym.Ones())*arg)
diff --git a/test/test_grudge.py b/test/test_grudge.py
index bc871da57335fd83a5bcac9d6e2cf4868cc1d2b2..3754f49f4b596b335004ab3c4574a9e2619b302c 100644
--- a/test/test_grudge.py
+++ b/test/test_grudge.py
@@ -85,6 +85,82 @@ def test_inverse_metric(ctx_getter, dim):
assert err < 1e-12, (i, j, err)
+def test_1d_mass_mat_trig(ctx_getter):
+ """Check the integral of some trig functions on an interval using the mass
+ matrix
+ """
+
+ cl_ctx = cl.create_some_context()
+ queue = cl.CommandQueue(cl_ctx)
+
+ from meshmode.mesh.generation import generate_regular_rect_mesh
+
+ mesh = generate_regular_rect_mesh(a=(-4*np.pi,), b=(9*np.pi,),
+ n=(17,), order=1)
+
+ discr = Discretization(cl_ctx, mesh, order=8)
+
+ x = sym.nodes(1)
+ f = bind(discr, sym.cos(x[0])**2)(queue)
+ ones = bind(discr, sym.Ones())(queue)
+
+ mass_op = bind(discr, sym.MassOperator()(sym.Field("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()
+
+ true_integral = 13*np.pi/2
+ err_1 = abs(num_integral_1-true_integral)
+ err_2 = abs(num_integral_2-true_integral)
+ err_3 = abs(num_integral_3-true_integral)
+
+ assert err_1 < 1e-10
+ assert err_2 < 1e-10
+ assert err_3 < 1e-10
+
+
+@pytest.mark.parametrize("dim", [1, 2, 3])
+def test_tri_diff_mat(ctx_getter, dim, order=4):
+ """Check differentiation matrix along the coordinate axes on a disk
+
+ Uses sines as the function to differentiate.
+ """
+
+ cl_ctx = cl.create_some_context()
+ queue = cl.CommandQueue(cl_ctx)
+
+ from meshmode.mesh.generation import generate_regular_rect_mesh
+
+ from pytools.convergence import EOCRecorder
+ axis_eoc_recs = [EOCRecorder() for axis in range(dim)]
+
+ for n in [10, 20]:
+ mesh = generate_regular_rect_mesh(a=(-0.5,)*dim, b=(0.5,)*dim,
+ n=(n,)*dim, order=4)
+
+ discr = Discretization(cl_ctx, mesh, order=4)
+ nabla = sym.nabla(dim)
+
+ for axis in range(dim):
+ x = sym.nodes(dim)
+
+ 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"))
+ bound_op = bind(discr, sym_op)
+ df_num = bound_op(queue, f=f)
+
+ linf_error = la.norm((df_num-df).get(), np.Inf)
+ axis_eoc_recs[axis].add_data_point(1/n, linf_error)
+
+ for axis, eoc_rec in enumerate(axis_eoc_recs):
+ print(axis)
+ print(eoc_rec)
+ assert eoc_rec.order_estimate() >= order
+
+
# You can test individual routines by typing
# $ python test_layer_pot.py 'test_routine()'