From f3d283606c6d85370110db1a9d6795454f5e8680 Mon Sep 17 00:00:00 2001 From: Alex Fikl Date: Tue, 19 Jun 2018 14:33:26 -0500 Subject: [PATCH 01/34] matrix: add a p2p matrix builder --- pytential/symbolic/matrix.py | 124 ++++++++++++++++++++++++++++------- 1 file changed, 100 insertions(+), 24 deletions(-) diff --git a/pytential/symbolic/matrix.py b/pytential/symbolic/matrix.py index f0b1d82f..c9d91859 100644 --- a/pytential/symbolic/matrix.py +++ b/pytential/symbolic/matrix.py @@ -34,23 +34,68 @@ import pytential.symbolic.primitives as sym from pytential.symbolic.execution import bind +# {{{ helpers + def is_zero(x): return isinstance(x, (int, float, complex, np.number)) and x == 0 +def _resample_arg(queue, source, x): + if not isinstance(x, np.ndarray): + return x + + if len(x.shape) >= 2: + raise RuntimeError("matrix variables in kernel arguments") + + def resample(y): + return source.resampler(queue, cl.array.to_device(queue, y)).get(queue) + + from pytools.obj_array import with_object_array_or_scalar + return with_object_array_or_scalar(resample, x) + + +def _get_layer_potential_args(mapper, expr, source): + kernel_args = {} + for arg_name, arg_expr in six.iteritems(expr.kernel_arguments): + rec_arg = mapper.rec(arg_expr) + kernel_args[arg_name] = _resample_arg(mapper.queue, source, rec_arg) + + return kernel_args + + +def _get_kernel_args(mapper, kernel, expr, source): + # NOTE: copied from pytential.symbolic.primitives.IntG + inner_kernel_args = kernel.get_args() + kernel.get_source_args() + inner_kernel_args = set(arg.loopy_arg.name for arg in inner_kernel_args) + + kernel_args = {} + for arg_name, arg_expr in six.iteritems(expr.kernel_arguments): + if arg_name not in inner_kernel_args: + continue + + rec_arg = mapper.rec(arg_expr) + kernel_args[arg_name] = _resample_arg(mapper.queue, source, rec_arg) + + return kernel_args + +# }}} + +# {{{ QBX layer potential matrix + # FIXME: PyOpenCL doesn't do all the required matrix math yet. # We'll cheat and build the matrix on the host. class MatrixBuilder(EvaluationMapperBase): def __init__(self, queue, dep_expr, other_dep_exprs, dep_source, places, context): + super(MatrixBuilder, self).__init__(context=context) + self.queue = queue self.dep_expr = dep_expr self.other_dep_exprs = other_dep_exprs self.dep_source = dep_source self.dep_discr = dep_source.density_discr self.places = places - self.context = context def map_variable(self, expr): if expr == self.dep_expr: @@ -152,27 +197,7 @@ class MatrixBuilder(EvaluationMapperBase): raise NotImplementedError("layer potentials on non-variables") kernel = expr.kernel - - kernel_args = {} - for arg_name, arg_expr in six.iteritems(expr.kernel_arguments): - rec_arg = self.rec(arg_expr) - - if isinstance(rec_arg, np.ndarray): - if len(rec_arg.shape) == 2: - raise RuntimeError("matrix variables in kernel arguments") - if len(rec_arg.shape) == 1: - from pytools.obj_array import with_object_array_or_scalar - - def resample(x): - return ( - source.resampler( - self.queue, - cl.array.to_device(self.queue, x)) - .get(queue=self.queue)) - - rec_arg = with_object_array_or_scalar(resample, rec_arg) - - kernel_args[arg_name] = rec_arg + kernel_args = _get_layer_potential_args(self, expr, source) from sumpy.expansion.local import LineTaylorLocalExpansion local_expn = LineTaylorLocalExpansion(kernel, source.qbx_order) @@ -181,8 +206,6 @@ class MatrixBuilder(EvaluationMapperBase): mat_gen = LayerPotentialMatrixGenerator( self.queue.context, (local_expn,)) - assert target_discr is source.density_discr - from pytential.qbx.utils import get_centers_on_side assert abs(expr.qbx_forced_limit) > 0 @@ -244,3 +267,56 @@ class MatrixBuilder(EvaluationMapperBase): result = result.get() return result + +# }}} + + +# {{{ p2p matrix + +class P2PMatrixBuilder(MatrixBuilder): + def __init__(self, queue, dep_expr, other_dep_exprs, dep_source, places, + context, exclude_self=True): + super(P2PMatrixBuilder, self).__init__(queue, + dep_expr, other_dep_exprs, dep_source, places, context) + + self.exclude_self = exclude_self + + def map_int_g(self, expr): + source = self.places[expr.source] + target_discr = self.places[expr.target] + + if source.density_discr is not target_discr: + raise NotImplementedError() + + rec_density = self.rec(expr.density) + if is_zero(rec_density): + return 0 + + assert isinstance(rec_density, np.ndarray) + if len(rec_density.shape) != 2: + raise NotImplementedError("layer potentials on non-variables") + + try: + kernel = expr.kernel.inner_kernel + except AttributeError: + kernel = expr.kernel + kernel_args = _get_kernel_args(self, kernel, expr, source) + if self.exclude_self: + kernel_args["target_to_source"] = np.arange(0, target_discr.nnodes) + + from sumpy.p2p import P2PMatrixGenerator + mat_gen = P2PMatrixGenerator( + self.queue.context, (kernel,), exclude_self=self.exclude_self) + + _, (mat,) = mat_gen(self.queue, + targets=target_discr.nodes(), + sources=source.density_discr.nodes(), + **kernel_args) + + mat = mat.get() + mat = mat.dot(rec_density) + + return mat +# }}} + +# vim: foldmethod=marker -- GitLab From b6ee542d0e6fa652991be353583fee374032ca5f Mon Sep 17 00:00:00 2001 From: Alex Fikl Date: Tue, 19 Jun 2018 14:56:44 -0500 Subject: [PATCH 02/34] matrix: add block matrix builders --- pytential/symbolic/matrix.py | 253 +++++++++++++++++++++++++++++++++-- 1 file changed, 245 insertions(+), 8 deletions(-) diff --git a/pytential/symbolic/matrix.py b/pytential/symbolic/matrix.py index c9d91859..5a62aa04 100644 --- a/pytential/symbolic/matrix.py +++ b/pytential/symbolic/matrix.py @@ -29,6 +29,8 @@ import pyopencl.array # noqa import six from six.moves import intern +from pytools import memoize_method + from pytential.symbolic.mappers import EvaluationMapperBase import pytential.symbolic.primitives as sym from pytential.symbolic.execution import bind @@ -41,6 +43,9 @@ def is_zero(x): def _resample_arg(queue, source, x): + if source is None: + return x + if not isinstance(x, np.ndarray): return x @@ -80,7 +85,8 @@ def _get_kernel_args(mapper, kernel, expr, source): # }}} -# {{{ QBX layer potential matrix + +# {{{ QBX layer potential matrix builder # FIXME: PyOpenCL doesn't do all the required matrix math yet. # We'll cheat and build the matrix on the host. @@ -207,15 +213,14 @@ class MatrixBuilder(EvaluationMapperBase): self.queue.context, (local_expn,)) from pytential.qbx.utils import get_centers_on_side - assert abs(expr.qbx_forced_limit) > 0 + _, (mat,) = mat_gen(self.queue, target_discr.nodes(), source.quad_stage2_density_discr.nodes(), get_centers_on_side(source, expr.qbx_forced_limit), expansion_radii=self.dep_source._expansion_radii("nsources"), **kernel_args) - mat = mat.get() waa = source.weights_and_area_elements().get(queue=self.queue) @@ -252,9 +257,7 @@ class MatrixBuilder(EvaluationMapperBase): arg, = expr.parameters rec_arg = self.rec(arg) - if ( - isinstance(rec_arg, np.ndarray) - and len(rec_arg.shape) == 2): + if isinstance(rec_arg, np.ndarray) and len(rec_arg.shape) == 2: raise RuntimeError("expression is nonlinear in variable") if isinstance(rec_arg, np.ndarray): @@ -271,7 +274,7 @@ class MatrixBuilder(EvaluationMapperBase): # }}} -# {{{ p2p matrix +# {{{ p2p matrix builder class P2PMatrixBuilder(MatrixBuilder): def __init__(self, queue, dep_expr, other_dep_exprs, dep_source, places, @@ -302,7 +305,8 @@ class P2PMatrixBuilder(MatrixBuilder): kernel = expr.kernel kernel_args = _get_kernel_args(self, kernel, expr, source) if self.exclude_self: - kernel_args["target_to_source"] = np.arange(0, target_discr.nnodes) + kernel_args["target_to_source"] = \ + cl.array.arange(self.queue, 0, target_discr.nnodes, dtype=np.int) from sumpy.p2p import P2PMatrixGenerator mat_gen = P2PMatrixGenerator( @@ -319,4 +323,237 @@ class P2PMatrixBuilder(MatrixBuilder): return mat # }}} + +# {{{ block matrix builders + +class MatrixBlockBuilderBase(EvaluationMapperBase): + def __init__(self, queue, dep_expr, other_dep_exprs, dep_source, + places, context, index_set): + super(MatrixBlockBuilderBase, self).__init__(context=context) + + self.queue = queue + self.dep_expr = dep_expr + self.other_dep_exprs = other_dep_exprs + self.dep_source = dep_source + self.dep_discr = dep_source.density_discr + self.places = places + + self.index_set = index_set + + def _map_dep_variable(self): + return np.eye(self.index_set.srcindices.shape[0]) + + def map_variable(self, expr): + if expr == self.dep_expr: + return self._map_dep_variable() + elif expr in self.other_dep_exprs: + return 0 + else: + return super(MatrixBlockBuilderBase, self).map_variable(expr) + + def map_subscript(self, expr): + if expr == self.dep_expr: + return self.variable_identity() + elif expr in self.other_dep_exprs: + return 0 + else: + return super(MatrixBlockBuilderBase, self).map_subscript(expr) + + def map_num_reference_derivative(self, expr): + rec_operand = self.rec(expr.operand) + + assert isinstance(rec_operand, np.ndarray) + if len(rec_operand.shape) == 2: + raise NotImplementedError("derivatives") + + where_discr = self.places[expr.where] + op = sym.NumReferenceDerivative(expr.ref_axes, sym.var("u")) + return bind(where_discr, op)( + self.queue, u=cl.array.to_device(self.queue, rec_operand)).get() + + def map_node_coordinate_component(self, expr): + where_discr = self.places[expr.where] + op = sym.NodeCoordinateComponent(expr.ambient_axis) + return bind(where_discr, op)(self.queue).get() + + def map_call(self, expr): + arg, = expr.parameters + rec_arg = self.rec(arg) + + if isinstance(rec_arg, np.ndarray) and len(rec_arg.shape) == 2: + raise RuntimeError("expression is nonlinear in variable") + + if isinstance(rec_arg, np.ndarray): + rec_arg = cl.array.to_device(self.queue, rec_arg) + + op = expr.function(sym.var("u")) + result = bind(self.dep_source, op)(self.queue, u=rec_arg) + + if isinstance(result, cl.array.Array): + result = result.get() + + return result + + +class MatrixBlockBuilder(MatrixBlockBuilderBase): + def __init__(self, queue, dep_expr, other_dep_exprs, dep_source, + places, context, index_set): + super(MatrixBlockBuilder, self).__init__(queue, + dep_expr, other_dep_exprs, dep_source, places, context, index_set) + + self.dummy = MatrixBlockBuilderBase(queue, + dep_expr, other_dep_exprs, dep_source, places, context, index_set) + + @memoize_method + def _oversampled_index_set(self, resampler): + from pytential.direct_solver import partition_from_coarse + srcindices, srcranges = partition_from_coarse(self.queue, resampler, + self.index_set.srcindices, self.index_set.srcranges) + + from sumpy.tools import MatrixBlockIndex + ovsm_index_set = MatrixBlockIndex(self.queue, + self.index_set.tgtindices, srcindices, + self.index_set.tgtranges, srcranges) + + return ovsm_index_set + + def _map_dep_variable(self): + return np.equal(self.index_set.tgtindices.reshape(-1, 1), + self.index_set.srcindices.reshape(1, -1)).astype(np.float64) + + def map_int_g(self, expr): + source = self.places[expr.source] + target_discr = self.places[expr.target] + + if source.density_discr is not target_discr: + raise NotImplementedError() + + rec_density = self.dummy.rec(expr.density) + if is_zero(rec_density): + return 0 + + assert isinstance(rec_density, np.ndarray) + if len(rec_density.shape) != 2: + raise NotImplementedError("layer potentials on non-variables") + + resampler = source.direct_resampler + ovsm_index_set = self._oversampled_index_set(resampler) + + kernel = expr.kernel + kernel_args = _get_layer_potential_args(self, expr, source) + + from sumpy.expansion.local import LineTaylorLocalExpansion + local_expn = LineTaylorLocalExpansion(kernel, source.qbx_order) + + from sumpy.qbx import LayerPotentialMatrixBlockGenerator + mat_gen = LayerPotentialMatrixBlockGenerator( + self.queue.context, (local_expn,)) + + from pytential.qbx.utils import get_centers_on_side + assert abs(expr.qbx_forced_limit) > 0 + + _, (mat,) = mat_gen(self.queue, + target_discr.nodes(), + source.quad_stage2_density_discr.nodes(), + get_centers_on_side(source, expr.qbx_forced_limit), + expansion_radii=self.dep_source._expansion_radii("nsources"), + index_set=ovsm_index_set, + **kernel_args) + mat = mat.get() + + ovsm_blkranges = ovsm_index_set.linear_ranges() + ovsm_rowindices, ovsm_colindices = ovsm_index_set.linear_indices() + + waa = source.weights_and_area_elements().get(queue=self.queue) + mat *= waa[ovsm_colindices] + + # TODO: there should be a better way to do the matmat required for + # the resampling + resample_mat = resampler.full_resample_matrix(self.queue).get(self.queue) + + rmat = np.empty(ovsm_blkranges.shape[0] - 1, dtype=np.object) + for i in range(ovsm_blkranges.shape[0] - 1): + # TODO: would be nice to move some of this to sumpy.MatrixBlockIndex + ovsm_iblk = np.s_[ovsm_blkranges[i]:ovsm_blkranges[i + 1]] + ovsm_shape = (ovsm_index_set.tgtranges[i + 1] - + ovsm_index_set.tgtranges[i], + ovsm_index_set.srcranges[i + 1] - + ovsm_index_set.srcranges[i]) + mat_blk = mat[ovsm_iblk].reshape(*ovsm_shape) + + itgt = np.s_[ovsm_index_set.srcranges[i]:ovsm_index_set.srcranges[i + 1]] + itgt = ovsm_index_set.srcindices[itgt] + isrc = np.s_[self.index_set.srcranges[i]:self.index_set.srcranges[i + 1]] + isrc = self.index_set.srcindices[isrc] + resample_mat_blk = resample_mat[np.ix_(itgt, isrc)] + + rmat[i] = mat_blk.dot(resample_mat_blk).reshape(-1) + mat = np.hstack(rmat) + + # TODO: multiply with rec_density + + return mat + + +class P2PMatrixBlockBuilder(MatrixBlockBuilderBase): + def __init__(self, queue, dep_expr, other_dep_exprs, dep_source, + places, context, index_set, exclude_self=True): + super(P2PMatrixBlockBuilder, self).__init__(queue, + dep_expr, other_dep_exprs, dep_source, places, context, index_set) + + self.dummy = MatrixBlockBuilderBase(queue, + dep_expr, other_dep_exprs, dep_source, places, context, index_set) + self.exclude_self = exclude_self + + def _map_dep_variable(self): + return np.equal(self.index_set.tgtindices.reshape(-1, 1), + self.index_set.srcindices.reshape(1, -1)).astype(np.float64) + + def block(self, i): + itgt = np.s_[self.index_set.tgtranges[i]:self.index_set.tgtranges[i + 1]] + isrc = np.s_[self.index_set.srcranges[i]:self.index_set.srcranges[i + 1]] + + return (itgt, isrc) + + def map_int_g(self, expr): + source = self.places[expr.source] + target_discr = self.places[expr.target] + + if source.density_discr is not target_discr: + raise NotImplementedError() + + rec_density = self.dummy.rec(expr.density) + if is_zero(rec_density): + return 0 + + assert isinstance(rec_density, np.ndarray) + if len(rec_density.shape) != 2: + raise NotImplementedError("layer potentials on non-variables") + + try: + kernel = expr.kernel.inner_kernel + except AttributeError: + kernel = expr.kernel + kernel_args = _get_kernel_args(self, kernel, expr, source) + if self.exclude_self: + kernel_args["target_to_source"] = \ + cl.array.arange(self.queue, 0, target_discr.nnodes, dtype=np.int) + + from sumpy.p2p import P2PMatrixBlockGenerator + mat_gen = P2PMatrixBlockGenerator( + self.queue.context, (kernel,), exclude_self=self.exclude_self) + + _, (mat,) = mat_gen(self.queue, + targets=target_discr.nodes(), + sources=source.nodes(), + index_set=self.index_set, + **kernel_args) + mat = mat.get() + + # TODO: need to multiply by rec_density + + return mat + +# }}} + # vim: foldmethod=marker -- GitLab From bed93c255556bce25ce7fef547b5dd01339280a0 Mon Sep 17 00:00:00 2001 From: Alex Fikl Date: Tue, 19 Jun 2018 15:21:36 -0500 Subject: [PATCH 03/34] matrix: add tests for the block matrix builders --- pytential/symbolic/matrix.py | 4 +- test/test_linalg_proxy.py | 1 + test/test_matrix.py | 180 ++++++++++++++++++++++++++++++++++- 3 files changed, 180 insertions(+), 5 deletions(-) diff --git a/pytential/symbolic/matrix.py b/pytential/symbolic/matrix.py index 5a62aa04..06e963ac 100644 --- a/pytential/symbolic/matrix.py +++ b/pytential/symbolic/matrix.py @@ -488,7 +488,7 @@ class MatrixBlockBuilder(MatrixBlockBuilderBase): resample_mat_blk = resample_mat[np.ix_(itgt, isrc)] rmat[i] = mat_blk.dot(resample_mat_blk).reshape(-1) - mat = np.hstack(rmat) + mat = np.hstack(rmat) # TODO: multiply with rec_density @@ -545,7 +545,7 @@ class P2PMatrixBlockBuilder(MatrixBlockBuilderBase): _, (mat,) = mat_gen(self.queue, targets=target_discr.nodes(), - sources=source.nodes(), + sources=source.density_discr.nodes(), index_set=self.index_set, **kernel_args) mat = mat.get() diff --git a/test/test_linalg_proxy.py b/test/test_linalg_proxy.py index e8a063ca..f40f4053 100644 --- a/test/test_linalg_proxy.py +++ b/test/test_linalg_proxy.py @@ -88,6 +88,7 @@ def _build_block_index(discr, if method == 'elements': factor = 1.0 + density_discr = qbx.density_discr if method == 'nodes': nnodes = discr.nnodes else: diff --git a/test/test_matrix.py b/test/test_matrix.py index 9f0ff0b4..cadb3d65 100644 --- a/test/test_matrix.py +++ b/test/test_matrix.py @@ -22,16 +22,21 @@ OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN THE SOFTWARE. """ +from functools import partial + import numpy as np import numpy.linalg as la + import pyopencl as cl -import pytest + +from sumpy.symbolic import USE_SYMENGINE from meshmode.mesh.generation import \ ellipse, NArmedStarfish, make_curve_mesh + from pytential import bind, sym -from functools import partial -from sumpy.symbolic import USE_SYMENGINE +from pytential.symbolic.primitives import DEFAULT_TARGET, DEFAULT_SOURCE +import pytest from pyopencl.tools import ( # noqa pytest_generate_tests_for_pyopencl as pytest_generate_tests) @@ -155,6 +160,175 @@ def test_matrix_build(ctx_factory, k, curve_f, layer_pot_id, visualize=False): assert rel_err < 1e-13 +@pytest.mark.parametrize("ndim", [2, 3]) +@pytest.mark.parametrize("factor", [1.0, 0.6]) +def test_p2p_block_builder(ctx_factory, factor, ndim, visualize=False): + ctx = ctx_factory() + queue = cl.CommandQueue(ctx) + + # prevent cache explosion + from sympy.core.cache import clear_cache + clear_cache() + + from test_direct_solver import _create_data, _create_indices + target_order = 2 if ndim == 3 else 7 + qbx, op, u_sym = _create_data(queue, target_order=target_order, ndim=ndim) + + nblks = 10 + srcindices, srcranges = _create_indices(qbx, nblks, + method='nodes', factor=factor, random=True) + tgtindices, tgtranges = _create_indices(qbx, nblks, + method='nodes', factor=factor, random=True) + nblks = srcranges.shape[0] - 1 + assert srcranges.shape[0] == tgtranges.shape[0] + + from pytential.symbolic.execution import prepare_places, prepare_expr + places = prepare_places(qbx) + expr = prepare_expr(places, op) + + from sumpy.tools import MatrixBlockIndex + index_set = MatrixBlockIndex(queue, + tgtindices, srcindices, tgtranges, srcranges) + + from pytential.symbolic.matrix import P2PMatrixBlockBuilder + mbuilder = P2PMatrixBlockBuilder(queue, + dep_expr=u_sym, + other_dep_exprs=[], + dep_source=places[DEFAULT_SOURCE], + places=places, + context={}, + index_set=index_set) + blk = mbuilder(expr) + + from pytential.symbolic.matrix import P2PMatrixBuilder + mbuilder = P2PMatrixBuilder(queue, + dep_expr=u_sym, + other_dep_exprs=[], + dep_source=places[DEFAULT_SOURCE], + places=places, + context={}) + mat = mbuilder(expr) + + if visualize and ndim == 2: + blk_full = np.zeros_like(mat) + mat_full = np.zeros_like(mat) + + blkranges = index_set.linear_ranges() + rowindices, colindices = index_set.linear_indices() + for i in range(srcranges.shape[0] - 1): + iblk = np.s_[blkranges[i]:blkranges[i + 1]] + itgt = rowindices[iblk] + isrc = colindices[iblk] + + blk_full[itgt, isrc] = blk[iblk] + mat_full[itgt, isrc] = mat[itgt, isrc] + + import matplotlib.pyplot as mp + _, (ax1, ax2) = mp.subplots(1, 2, + figsize=(10, 8), dpi=300, constrained_layout=True) + ax1.imshow(blk_full) + ax1.set_title('P2PMatrixBlockBuilder') + ax2.imshow(mat_full) + ax2.set_title('P2PMatrixBuilder') + mp.savefig("test_p2p_block_{}d_{:.1f}.png".format(ndim, factor)) + + blkranges = index_set.linear_ranges() + rowindices, colindices = index_set.linear_indices() + for i in range(nblks): + iblk = np.s_[blkranges[i]:blkranges[i + 1]] + itgt = rowindices[iblk] + isrc = colindices[iblk] + + eps = 1.0e-14 * la.norm(mat[itgt, isrc]) + if visualize: + print('block[{:04}]: {:.5e}'.format(i, + la.norm(blk[iblk] - mat[itgt, isrc].reshape(-1)))) + assert la.norm(blk[iblk] - mat[itgt, isrc].reshape(-1)) < eps + + +@pytest.mark.parametrize("ndim", [2, 3]) +def test_qbx_block_builder(ctx_factory, ndim, visualize=False): + ctx = ctx_factory() + queue = cl.CommandQueue(ctx) + + # prevent cache explosion + from sympy.core.cache import clear_cache + clear_cache() + + from test_direct_solver import _create_data, _create_indices + target_order = 2 if ndim == 3 else 7 + qbx, op, u_sym = _create_data(queue, target_order=target_order, ndim=ndim) + + nblks = 10 + tgtindices, tgtranges = _create_indices(qbx, nblks) + srcindices, srcranges = _create_indices(qbx, nblks) + nblks = srcranges.shape[0] - 1 + assert srcranges.shape[0] == tgtranges.shape[0] + + from pytential.symbolic.execution import prepare_places, prepare_expr + places = prepare_places(qbx) + expr = prepare_expr(places, op) + + from sumpy.tools import MatrixBlockIndex + index_set = MatrixBlockIndex(queue, + tgtindices, srcindices, tgtranges, srcranges) + + from pytential.symbolic.matrix import MatrixBlockBuilder + mbuilder = MatrixBlockBuilder(queue, + dep_expr=u_sym, + other_dep_exprs=[], + dep_source=places[DEFAULT_SOURCE], + places=places, + context={}, + index_set=index_set) + blk = mbuilder(expr) + + from pytential.symbolic.matrix import MatrixBuilder + mbuilder = MatrixBuilder(queue, + dep_expr=u_sym, + other_dep_exprs=[], + dep_source=places[DEFAULT_SOURCE], + places=places, + context={}) + mat = mbuilder(expr) + + if visualize: + blk_full = np.zeros_like(mat) + mat_full = np.zeros_like(mat) + + blkranges = index_set.linear_ranges() + rowindices, colindices = index_set.linear_indices() + for i in range(srcranges.shape[0] - 1): + iblk = np.s_[blkranges[i]:blkranges[i + 1]] + itgt = rowindices[iblk] + isrc = colindices[iblk] + + blk_full[itgt, isrc] = blk[iblk] + mat_full[itgt, isrc] = mat[itgt, isrc] + + import matplotlib.pyplot as mp + _, (ax1, ax2) = mp.subplots(1, 2, + figsize=(10, 8), constrained_layout=True) + ax1.imshow(mat_full) + ax1.set_title('MatrixBuilder') + ax2.imshow(blk_full) + ax2.set_title('MatrixBlockBuilder') + mp.savefig("test_qbx_block_builder.png", dpi=300) + + blkranges = index_set.linear_ranges() + rowindices, colindices = index_set.linear_indices() + for i in range(nblks): + iblk = np.s_[blkranges[i]:blkranges[i + 1]] + itgt = rowindices[iblk] + isrc = colindices[iblk] + + eps = 1.0e-14 * la.norm(mat[itgt, isrc]) + if visualize: + print('block[{:04}]: {:.5e}'.format(i, + la.norm(blk[iblk] - mat[itgt, isrc].reshape(-1)))) + assert la.norm(blk[iblk] - mat[itgt, isrc].reshape(-1)) < eps + + if __name__ == "__main__": import sys if len(sys.argv) > 1: -- GitLab From 0bc2500a12da07dc7822a919f305e1a485a51824 Mon Sep 17 00:00:00 2001 From: Alex Fikl Date: Tue, 19 Jun 2018 15:25:22 -0500 Subject: [PATCH 04/34] tests: flake8 --- test/test_matrix.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/test/test_matrix.py b/test/test_matrix.py index cadb3d65..ed440632 100644 --- a/test/test_matrix.py +++ b/test/test_matrix.py @@ -34,7 +34,7 @@ from meshmode.mesh.generation import \ ellipse, NArmedStarfish, make_curve_mesh from pytential import bind, sym -from pytential.symbolic.primitives import DEFAULT_TARGET, DEFAULT_SOURCE +from pytential.symbolic.primitives import DEFAULT_SOURCE import pytest from pyopencl.tools import ( # noqa -- GitLab From 3aefe76d66df6daadf5ebb06ca754058ffe7ceef Mon Sep 17 00:00:00 2001 From: Alex Fikl Date: Sun, 24 Jun 2018 19:49:28 -0500 Subject: [PATCH 05/34] direct-solver: modify to latest changes in linalg/proxy --- pytential/symbolic/matrix.py | 9 ++- test/test_linalg_proxy.py | 1 - test/test_matrix.py | 119 ++++++++++++++++++++--------------- 3 files changed, 75 insertions(+), 54 deletions(-) diff --git a/pytential/symbolic/matrix.py b/pytential/symbolic/matrix.py index 06e963ac..034b5ac8 100644 --- a/pytential/symbolic/matrix.py +++ b/pytential/symbolic/matrix.py @@ -406,14 +406,17 @@ class MatrixBlockBuilder(MatrixBlockBuilderBase): @memoize_method def _oversampled_index_set(self, resampler): - from pytential.direct_solver import partition_from_coarse + from pytential.linalg.proxy import partition_from_coarse srcindices, srcranges = partition_from_coarse(self.queue, resampler, self.index_set.srcindices, self.index_set.srcranges) from sumpy.tools import MatrixBlockIndex + # TODO: fix MatrixBlockIndex to work with CL arrays ovsm_index_set = MatrixBlockIndex(self.queue, - self.index_set.tgtindices, srcindices, - self.index_set.tgtranges, srcranges) + self.index_set.tgtindices, + srcindices.get(self.queue), + self.index_set.tgtranges, + srcranges.get(self.queue)) return ovsm_index_set diff --git a/test/test_linalg_proxy.py b/test/test_linalg_proxy.py index f40f4053..e8a063ca 100644 --- a/test/test_linalg_proxy.py +++ b/test/test_linalg_proxy.py @@ -88,7 +88,6 @@ def _build_block_index(discr, if method == 'elements': factor = 1.0 - density_discr = qbx.density_discr if method == 'nodes': nnodes = discr.nnodes else: diff --git a/test/test_matrix.py b/test/test_matrix.py index ed440632..b6810730 100644 --- a/test/test_matrix.py +++ b/test/test_matrix.py @@ -42,38 +42,25 @@ from pyopencl.tools import ( # noqa as pytest_generate_tests) -@pytest.mark.skipif(USE_SYMENGINE, - reason="https://gitlab.tiker.net/inducer/sumpy/issues/25") -@pytest.mark.parametrize("k", [0, 42]) -@pytest.mark.parametrize("curve_f", [ - partial(ellipse, 3), - NArmedStarfish(5, 0.25)]) -@pytest.mark.parametrize("layer_pot_id", [1, 2]) -def test_matrix_build(ctx_factory, k, curve_f, layer_pot_id, visualize=False): - cl_ctx = ctx_factory() - queue = cl.CommandQueue(cl_ctx) - - # prevent cache 'splosion - from sympy.core.cache import clear_cache - clear_cache() - - target_order = 7 - qbx_order = 4 - nelements = 30 - +def _build_op(lpot_id, k=0, ndim=2): from sumpy.kernel import LaplaceKernel, HelmholtzKernel + knl_kwargs = {"qbx_forced_limit": "avg"} if k: - knl = HelmholtzKernel(2) - knl_kwargs = {"k": k} + knl = HelmholtzKernel(ndim) + knl_kwargs["k"] = k else: - knl = LaplaceKernel(2) - knl_kwargs = {} + knl = LaplaceKernel(ndim) - from pytools.obj_array import make_obj_array, is_obj_array - if layer_pot_id == 1: + if lpot_id == 1: + # scalar single-layer potential + u_sym = sym.var("u") + op = sym.S(knl, u_sym, **knl_kwargs) + elif lpot_id == 2: + # scalar double-layer potential u_sym = sym.var("u") - op = sym.Sp(knl, u_sym, **knl_kwargs) - elif layer_pot_id == 2: + op = sym.D(knl, u_sym, **knl_kwargs) + elif lpot_id == 3: + # vector potential u_sym = sym.make_sym_vector("u", 2) u0_sym, u1_sym = u_sym @@ -85,8 +72,29 @@ def test_matrix_build(ctx_factory, k, curve_f, layer_pot_id, visualize=False): 0.3 * sym.D(knl, u0_sym, **knl_kwargs) ]) else: - raise ValueError("Unknown layer_pot_id: {}".format(layer_pot_id)) + raise ValueError("Unknown lpot_id: {}".format(lpot_id)) + + return op, u_sym + +@pytest.mark.skipif(USE_SYMENGINE, + reason="https://gitlab.tiker.net/inducer/sumpy/issues/25") +@pytest.mark.parametrize("k", [0, 42]) +@pytest.mark.parametrize("curve_f", [ + partial(ellipse, 3), + NArmedStarfish(5, 0.25)]) +@pytest.mark.parametrize("lpot_id", [1, 3]) +def test_matrix_build(ctx_factory, k, curve_f, lpot_id, visualize=False): + cl_ctx = ctx_factory() + queue = cl.CommandQueue(cl_ctx) + + # prevent cache 'splosion + from sympy.core.cache import clear_cache + clear_cache() + + target_order = 7 + qbx_order = 4 + nelements = 30 mesh = make_curve_mesh(curve_f, np.linspace(0, 1, nelements + 1), target_order) @@ -94,16 +102,18 @@ def test_matrix_build(ctx_factory, k, curve_f, layer_pot_id, visualize=False): from meshmode.discretization import Discretization from meshmode.discretization.poly_element import \ InterpolatoryQuadratureSimplexGroupFactory - from pytential.qbx import QBXLayerPotentialSource pre_density_discr = Discretization( cl_ctx, mesh, InterpolatoryQuadratureSimplexGroupFactory(target_order)) + from pytential.qbx import QBXLayerPotentialSource qbx, _ = QBXLayerPotentialSource(pre_density_discr, 4 * target_order, qbx_order, # Don't use FMM for now fmm_order=False).with_refinement() density_discr = qbx.density_discr + + op, u_sym = _build_op(lpot_id, k=k) bound_op = bind(qbx, op) from pytential.symbolic.execution import build_matrix @@ -162,7 +172,9 @@ def test_matrix_build(ctx_factory, k, curve_f, layer_pot_id, visualize=False): @pytest.mark.parametrize("ndim", [2, 3]) @pytest.mark.parametrize("factor", [1.0, 0.6]) -def test_p2p_block_builder(ctx_factory, factor, ndim, visualize=False): +@pytest.mark.parametrize("lpot_id", [1, 2]) +def test_p2p_block_builder(ctx_factory, factor, ndim, lpot_id, + visualize=False): ctx = ctx_factory() queue = cl.CommandQueue(ctx) @@ -170,16 +182,15 @@ def test_p2p_block_builder(ctx_factory, factor, ndim, visualize=False): from sympy.core.cache import clear_cache clear_cache() - from test_direct_solver import _create_data, _create_indices + from test_linalg_proxy import _build_qbx_discr, _build_block_index target_order = 2 if ndim == 3 else 7 - qbx, op, u_sym = _create_data(queue, target_order=target_order, ndim=ndim) - - nblks = 10 - srcindices, srcranges = _create_indices(qbx, nblks, - method='nodes', factor=factor, random=True) - tgtindices, tgtranges = _create_indices(qbx, nblks, - method='nodes', factor=factor, random=True) - nblks = srcranges.shape[0] - 1 + qbx = _build_qbx_discr(queue, target_order=target_order, ndim=ndim) + op, u_sym = _build_op(lpot_id, ndim=ndim) + + srcindices, srcranges = _build_block_index(queue, qbx.density_discr, + method='nodes', factor=factor) + tgtindices, tgtranges = _build_block_index(queue, qbx.density_discr, + method='nodes', factor=factor) assert srcranges.shape[0] == tgtranges.shape[0] from pytential.symbolic.execution import prepare_places, prepare_expr @@ -187,8 +198,12 @@ def test_p2p_block_builder(ctx_factory, factor, ndim, visualize=False): expr = prepare_expr(places, op) from sumpy.tools import MatrixBlockIndex + # TODO: make MatrixBlockIndex work properly with CL arrays index_set = MatrixBlockIndex(queue, - tgtindices, srcindices, tgtranges, srcranges) + tgtindices.get(queue), + srcindices.get(queue), + tgtranges.get(queue), + srcranges.get(queue)) from pytential.symbolic.matrix import P2PMatrixBlockBuilder mbuilder = P2PMatrixBlockBuilder(queue, @@ -234,7 +249,7 @@ def test_p2p_block_builder(ctx_factory, factor, ndim, visualize=False): blkranges = index_set.linear_ranges() rowindices, colindices = index_set.linear_indices() - for i in range(nblks): + for i in range(blkranges.shape[0] - 1): iblk = np.s_[blkranges[i]:blkranges[i + 1]] itgt = rowindices[iblk] isrc = colindices[iblk] @@ -247,7 +262,8 @@ def test_p2p_block_builder(ctx_factory, factor, ndim, visualize=False): @pytest.mark.parametrize("ndim", [2, 3]) -def test_qbx_block_builder(ctx_factory, ndim, visualize=False): +@pytest.mark.parametrize("lpot_id", [1]) +def test_qbx_block_builder(ctx_factory, ndim, lpot_id, visualize=False): ctx = ctx_factory() queue = cl.CommandQueue(ctx) @@ -255,14 +271,13 @@ def test_qbx_block_builder(ctx_factory, ndim, visualize=False): from sympy.core.cache import clear_cache clear_cache() - from test_direct_solver import _create_data, _create_indices + from test_linalg_proxy import _build_qbx_discr, _build_block_index target_order = 2 if ndim == 3 else 7 - qbx, op, u_sym = _create_data(queue, target_order=target_order, ndim=ndim) + qbx = _build_qbx_discr(queue, target_order=target_order, ndim=ndim) + op, u_sym = _build_op(lpot_id, ndim=ndim) - nblks = 10 - tgtindices, tgtranges = _create_indices(qbx, nblks) - srcindices, srcranges = _create_indices(qbx, nblks) - nblks = srcranges.shape[0] - 1 + tgtindices, tgtranges = _build_block_index(queue, qbx.density_discr) + srcindices, srcranges = _build_block_index(queue, qbx.density_discr) assert srcranges.shape[0] == tgtranges.shape[0] from pytential.symbolic.execution import prepare_places, prepare_expr @@ -270,8 +285,12 @@ def test_qbx_block_builder(ctx_factory, ndim, visualize=False): expr = prepare_expr(places, op) from sumpy.tools import MatrixBlockIndex + # TODO: make MatrixBlockIndex work properly with CL arrays index_set = MatrixBlockIndex(queue, - tgtindices, srcindices, tgtranges, srcranges) + tgtindices.get(queue), + srcindices.get(queue), + tgtranges.get(queue), + srcranges.get(queue)) from pytential.symbolic.matrix import MatrixBlockBuilder mbuilder = MatrixBlockBuilder(queue, @@ -317,7 +336,7 @@ def test_qbx_block_builder(ctx_factory, ndim, visualize=False): blkranges = index_set.linear_ranges() rowindices, colindices = index_set.linear_indices() - for i in range(nblks): + for i in range(blkranges.shape[0] - 1): iblk = np.s_[blkranges[i]:blkranges[i + 1]] itgt = rowindices[iblk] isrc = colindices[iblk] -- GitLab From 08e4441a6abec8a5b2f5eda852923ed01d7d0afd Mon Sep 17 00:00:00 2001 From: Alex Fikl Date: Sun, 24 Jun 2018 20:27:08 -0500 Subject: [PATCH 06/34] tests: add missing import --- test/test_matrix.py | 14 +++++++++----- 1 file changed, 9 insertions(+), 5 deletions(-) diff --git a/test/test_matrix.py b/test/test_matrix.py index b6810730..457de6ca 100644 --- a/test/test_matrix.py +++ b/test/test_matrix.py @@ -1,6 +1,9 @@ from __future__ import division, absolute_import, print_function -__copyright__ = "Copyright (C) 2015 Andreas Kloeckner" +__copyright__ = """ +Copyright (C) 2015 Andreas Kloeckner +Copyright (C) 2018 Andreas Kloeckner +""" __license__ = """ Permission is hereby granted, free of charge, to any person obtaining a copy @@ -28,6 +31,7 @@ import numpy as np import numpy.linalg as la import pyopencl as cl +from pytools.obj_array import make_obj_array, is_obj_array from sumpy.symbolic import USE_SYMENGINE from meshmode.mesh.generation import \ @@ -74,7 +78,7 @@ def _build_op(lpot_id, k=0, ndim=2): else: raise ValueError("Unknown lpot_id: {}".format(lpot_id)) - return op, u_sym + return op, u_sym, knl_kwargs @pytest.mark.skipif(USE_SYMENGINE, @@ -113,7 +117,7 @@ def test_matrix_build(ctx_factory, k, curve_f, lpot_id, visualize=False): fmm_order=False).with_refinement() density_discr = qbx.density_discr - op, u_sym = _build_op(lpot_id, k=k) + op, u_sym, knl_kwargs = _build_op(lpot_id, k=k) bound_op = bind(qbx, op) from pytential.symbolic.execution import build_matrix @@ -185,7 +189,7 @@ def test_p2p_block_builder(ctx_factory, factor, ndim, lpot_id, from test_linalg_proxy import _build_qbx_discr, _build_block_index target_order = 2 if ndim == 3 else 7 qbx = _build_qbx_discr(queue, target_order=target_order, ndim=ndim) - op, u_sym = _build_op(lpot_id, ndim=ndim) + op, u_sym, _ = _build_op(lpot_id, ndim=ndim) srcindices, srcranges = _build_block_index(queue, qbx.density_discr, method='nodes', factor=factor) @@ -274,7 +278,7 @@ def test_qbx_block_builder(ctx_factory, ndim, lpot_id, visualize=False): from test_linalg_proxy import _build_qbx_discr, _build_block_index target_order = 2 if ndim == 3 else 7 qbx = _build_qbx_discr(queue, target_order=target_order, ndim=ndim) - op, u_sym = _build_op(lpot_id, ndim=ndim) + op, u_sym, _ = _build_op(lpot_id, ndim=ndim) tgtindices, tgtranges = _build_block_index(queue, qbx.density_discr) srcindices, srcranges = _build_block_index(queue, qbx.density_discr) -- GitLab From 9dcdbed4a315f5f4fa1945af836b6dad665eb045 Mon Sep 17 00:00:00 2001 From: Alex Fikl Date: Fri, 29 Jun 2018 15:01:26 -0500 Subject: [PATCH 07/34] port to new modifications in MatrixBlockRanges --- pytential/symbolic/matrix.py | 70 +++++++++++----------------- test/test_matrix.py | 88 +++++++++++++----------------------- 2 files changed, 57 insertions(+), 101 deletions(-) diff --git a/pytential/symbolic/matrix.py b/pytential/symbolic/matrix.py index 034b5ac8..c534d64a 100644 --- a/pytential/symbolic/matrix.py +++ b/pytential/symbolic/matrix.py @@ -341,7 +341,7 @@ class MatrixBlockBuilderBase(EvaluationMapperBase): self.index_set = index_set def _map_dep_variable(self): - return np.eye(self.index_set.srcindices.shape[0]) + return np.eye(self.index_set.col.indices.shape[0]) def map_variable(self, expr): if expr == self.dep_expr: @@ -407,22 +407,19 @@ class MatrixBlockBuilder(MatrixBlockBuilderBase): @memoize_method def _oversampled_index_set(self, resampler): from pytential.linalg.proxy import partition_from_coarse - srcindices, srcranges = partition_from_coarse(self.queue, resampler, - self.index_set.srcindices, self.index_set.srcranges) + srcindices = partition_from_coarse(resampler, self.index_set) - from sumpy.tools import MatrixBlockIndex - # TODO: fix MatrixBlockIndex to work with CL arrays - ovsm_index_set = MatrixBlockIndex(self.queue, - self.index_set.tgtindices, - srcindices.get(self.queue), - self.index_set.tgtranges, - srcranges.get(self.queue)) + from sumpy.tools import MatrixBlockIndexRanges + ovsm_index_set = MatrixBlockIndexRanges(self.queue.context, + self.index_set.row, srcindices) return ovsm_index_set def _map_dep_variable(self): - return np.equal(self.index_set.tgtindices.reshape(-1, 1), - self.index_set.srcindices.reshape(1, -1)).astype(np.float64) + tgtindices = self.index_set.row.indices.get(self.queue).reshape(-1, 1) + srcindices = self.index_set.col.indices.get(self.queue).reshape(1, -1) + + return np.equal(tgtindices, srcindices).astype(np.float64) def map_int_g(self, expr): source = self.places[expr.source] @@ -462,38 +459,27 @@ class MatrixBlockBuilder(MatrixBlockBuilderBase): expansion_radii=self.dep_source._expansion_radii("nsources"), index_set=ovsm_index_set, **kernel_args) - mat = mat.get() - ovsm_blkranges = ovsm_index_set.linear_ranges() - ovsm_rowindices, ovsm_colindices = ovsm_index_set.linear_indices() + waa = source.weights_and_area_elements().with_queue(self.queue) + mat *= waa[ovsm_index_set.linear_col_indices] - waa = source.weights_and_area_elements().get(queue=self.queue) - mat *= waa[ovsm_colindices] + from sumpy.tools import MatrixBlockIndexRanges + ovsm_col_index = ovsm_index_set.get(self.queue) + ovsm_row_index = MatrixBlockIndexRanges(ovsm_col_index.cl_context, + ovsm_col_index.col, ovsm_col_index.row) - # TODO: there should be a better way to do the matmat required for - # the resampling + mat = mat.get(self.queue) resample_mat = resampler.full_resample_matrix(self.queue).get(self.queue) - rmat = np.empty(ovsm_blkranges.shape[0] - 1, dtype=np.object) - for i in range(ovsm_blkranges.shape[0] - 1): - # TODO: would be nice to move some of this to sumpy.MatrixBlockIndex - ovsm_iblk = np.s_[ovsm_blkranges[i]:ovsm_blkranges[i + 1]] - ovsm_shape = (ovsm_index_set.tgtranges[i + 1] - - ovsm_index_set.tgtranges[i], - ovsm_index_set.srcranges[i + 1] - - ovsm_index_set.srcranges[i]) - mat_blk = mat[ovsm_iblk].reshape(*ovsm_shape) - - itgt = np.s_[ovsm_index_set.srcranges[i]:ovsm_index_set.srcranges[i + 1]] - itgt = ovsm_index_set.srcindices[itgt] - isrc = np.s_[self.index_set.srcranges[i]:self.index_set.srcranges[i + 1]] - isrc = self.index_set.srcindices[isrc] - resample_mat_blk = resample_mat[np.ix_(itgt, isrc)] - - rmat[i] = mat_blk.dot(resample_mat_blk).reshape(-1) + rmat = np.empty(ovsm_index_set.nblocks, dtype=np.object) + for i in range(ovsm_index_set.nblocks): + mat_blk = ovsm_col_index.block_take(mat, i) + resample_blk = ovsm_row_index.take(resample_mat, i) + + rmat[i] = mat_blk.dot(resample_blk).reshape(-1) mat = np.hstack(rmat) - # TODO: multiply with rec_density + # TODO:: multiply with rec_density return mat @@ -509,14 +495,10 @@ class P2PMatrixBlockBuilder(MatrixBlockBuilderBase): self.exclude_self = exclude_self def _map_dep_variable(self): - return np.equal(self.index_set.tgtindices.reshape(-1, 1), - self.index_set.srcindices.reshape(1, -1)).astype(np.float64) - - def block(self, i): - itgt = np.s_[self.index_set.tgtranges[i]:self.index_set.tgtranges[i + 1]] - isrc = np.s_[self.index_set.srcranges[i]:self.index_set.srcranges[i + 1]] + tgtindices = self.index_set.row.indices.get(self.queue).reshape(-1, 1) + srcindices = self.index_set.col.indices.get(self.queue).reshape(1, -1) - return (itgt, isrc) + return np.equal(tgtindices, srcindices).astype(np.float64) def map_int_g(self, expr): source = self.places[expr.source] diff --git a/test/test_matrix.py b/test/test_matrix.py index 457de6ca..3becd048 100644 --- a/test/test_matrix.py +++ b/test/test_matrix.py @@ -2,7 +2,7 @@ from __future__ import division, absolute_import, print_function __copyright__ = """ Copyright (C) 2015 Andreas Kloeckner -Copyright (C) 2018 Andreas Kloeckner +Copyright (C) 2018 Alexandru Fikl """ __license__ = """ @@ -191,23 +191,17 @@ def test_p2p_block_builder(ctx_factory, factor, ndim, lpot_id, qbx = _build_qbx_discr(queue, target_order=target_order, ndim=ndim) op, u_sym, _ = _build_op(lpot_id, ndim=ndim) - srcindices, srcranges = _build_block_index(queue, qbx.density_discr, + srcindices = _build_block_index(qbx.density_discr, method='nodes', factor=factor) - tgtindices, tgtranges = _build_block_index(queue, qbx.density_discr, + tgtindices = _build_block_index(qbx.density_discr, method='nodes', factor=factor) - assert srcranges.shape[0] == tgtranges.shape[0] from pytential.symbolic.execution import prepare_places, prepare_expr places = prepare_places(qbx) expr = prepare_expr(places, op) - from sumpy.tools import MatrixBlockIndex - # TODO: make MatrixBlockIndex work properly with CL arrays - index_set = MatrixBlockIndex(queue, - tgtindices.get(queue), - srcindices.get(queue), - tgtranges.get(queue), - srcranges.get(queue)) + from sumpy.tools import MatrixBlockIndexRanges + index_set = MatrixBlockIndexRanges(ctx, tgtindices, srcindices) from pytential.symbolic.matrix import P2PMatrixBlockBuilder mbuilder = P2PMatrixBlockBuilder(queue, @@ -228,19 +222,16 @@ def test_p2p_block_builder(ctx_factory, factor, ndim, lpot_id, context={}) mat = mbuilder(expr) + index_set = index_set.get(queue) if visualize and ndim == 2: blk_full = np.zeros_like(mat) mat_full = np.zeros_like(mat) - blkranges = index_set.linear_ranges() - rowindices, colindices = index_set.linear_indices() - for i in range(srcranges.shape[0] - 1): - iblk = np.s_[blkranges[i]:blkranges[i + 1]] - itgt = rowindices[iblk] - isrc = colindices[iblk] + for i in range(index_set.nblocks): + itgt, isrc = index_set.block_indices(i) - blk_full[itgt, isrc] = blk[iblk] - mat_full[itgt, isrc] = mat[itgt, isrc] + blk_full[np.ix_(itgt, isrc)] = index_set.block_take(blk, i) + mat_full[np.ix_(itgt, isrc)] = index_set.take(mat, i) import matplotlib.pyplot as mp _, (ax1, ax2) = mp.subplots(1, 2, @@ -251,18 +242,14 @@ def test_p2p_block_builder(ctx_factory, factor, ndim, lpot_id, ax2.set_title('P2PMatrixBuilder') mp.savefig("test_p2p_block_{}d_{:.1f}.png".format(ndim, factor)) - blkranges = index_set.linear_ranges() - rowindices, colindices = index_set.linear_indices() - for i in range(blkranges.shape[0] - 1): - iblk = np.s_[blkranges[i]:blkranges[i + 1]] - itgt = rowindices[iblk] - isrc = colindices[iblk] + for i in range(index_set.nblocks): + eps = 1.0e-14 * la.norm(index_set.take(mat, i)) + error = la.norm(index_set.block_take(blk, i) - + index_set.take(mat, i)) - eps = 1.0e-14 * la.norm(mat[itgt, isrc]) if visualize: - print('block[{:04}]: {:.5e}'.format(i, - la.norm(blk[iblk] - mat[itgt, isrc].reshape(-1)))) - assert la.norm(blk[iblk] - mat[itgt, isrc].reshape(-1)) < eps + print('block[{:04}]: {:.5e}'.format(i, error)) + assert error < eps @pytest.mark.parametrize("ndim", [2, 3]) @@ -280,21 +267,15 @@ def test_qbx_block_builder(ctx_factory, ndim, lpot_id, visualize=False): qbx = _build_qbx_discr(queue, target_order=target_order, ndim=ndim) op, u_sym, _ = _build_op(lpot_id, ndim=ndim) - tgtindices, tgtranges = _build_block_index(queue, qbx.density_discr) - srcindices, srcranges = _build_block_index(queue, qbx.density_discr) - assert srcranges.shape[0] == tgtranges.shape[0] + tgtindices = _build_block_index(qbx.density_discr) + srcindices = _build_block_index(qbx.density_discr) from pytential.symbolic.execution import prepare_places, prepare_expr places = prepare_places(qbx) expr = prepare_expr(places, op) - from sumpy.tools import MatrixBlockIndex - # TODO: make MatrixBlockIndex work properly with CL arrays - index_set = MatrixBlockIndex(queue, - tgtindices.get(queue), - srcindices.get(queue), - tgtranges.get(queue), - srcranges.get(queue)) + from sumpy.tools import MatrixBlockIndexRanges + index_set = MatrixBlockIndexRanges(ctx, tgtindices, srcindices) from pytential.symbolic.matrix import MatrixBlockBuilder mbuilder = MatrixBlockBuilder(queue, @@ -315,19 +296,16 @@ def test_qbx_block_builder(ctx_factory, ndim, lpot_id, visualize=False): context={}) mat = mbuilder(expr) + index_set = index_set.get(queue) if visualize: blk_full = np.zeros_like(mat) mat_full = np.zeros_like(mat) - blkranges = index_set.linear_ranges() - rowindices, colindices = index_set.linear_indices() - for i in range(srcranges.shape[0] - 1): - iblk = np.s_[blkranges[i]:blkranges[i + 1]] - itgt = rowindices[iblk] - isrc = colindices[iblk] + for i in range(index_set.nblocks): + itgt, isrc = index_set.block_indices(i) - blk_full[itgt, isrc] = blk[iblk] - mat_full[itgt, isrc] = mat[itgt, isrc] + blk_full[np.ix_(itgt, isrc)] = index_set.block_take(blk, i) + mat_full[np.ix_(itgt, isrc)] = index_set.take(mat, i) import matplotlib.pyplot as mp _, (ax1, ax2) = mp.subplots(1, 2, @@ -338,18 +316,14 @@ def test_qbx_block_builder(ctx_factory, ndim, lpot_id, visualize=False): ax2.set_title('MatrixBlockBuilder') mp.savefig("test_qbx_block_builder.png", dpi=300) - blkranges = index_set.linear_ranges() - rowindices, colindices = index_set.linear_indices() - for i in range(blkranges.shape[0] - 1): - iblk = np.s_[blkranges[i]:blkranges[i + 1]] - itgt = rowindices[iblk] - isrc = colindices[iblk] + for i in range(index_set.nblocks): + eps = 1.0e-14 * la.norm(index_set.take(mat, i)) + error = la.norm(index_set.block_take(blk, i) - + index_set.take(mat, i)) - eps = 1.0e-14 * la.norm(mat[itgt, isrc]) if visualize: - print('block[{:04}]: {:.5e}'.format(i, - la.norm(blk[iblk] - mat[itgt, isrc].reshape(-1)))) - assert la.norm(blk[iblk] - mat[itgt, isrc].reshape(-1)) < eps + print('block[{:04}]: {:.5e}'.format(i, error)) + assert error < eps if __name__ == "__main__": -- GitLab From 766aec716bb6344a40d7bf1d3d668ef8a276bacc Mon Sep 17 00:00:00 2001 From: Alex Fikl Date: Tue, 3 Jul 2018 13:48:56 -0500 Subject: [PATCH 08/34] direct-solver: rename block builders --- pytential/symbolic/matrix.py | 8 ++++---- test/test_matrix.py | 12 ++++++------ 2 files changed, 10 insertions(+), 10 deletions(-) diff --git a/pytential/symbolic/matrix.py b/pytential/symbolic/matrix.py index c534d64a..4ea90d11 100644 --- a/pytential/symbolic/matrix.py +++ b/pytential/symbolic/matrix.py @@ -395,10 +395,10 @@ class MatrixBlockBuilderBase(EvaluationMapperBase): return result -class MatrixBlockBuilder(MatrixBlockBuilderBase): +class NearFieldBlockBuilder(MatrixBlockBuilderBase): def __init__(self, queue, dep_expr, other_dep_exprs, dep_source, places, context, index_set): - super(MatrixBlockBuilder, self).__init__(queue, + super(NearFieldBlockBuilder, self).__init__(queue, dep_expr, other_dep_exprs, dep_source, places, context, index_set) self.dummy = MatrixBlockBuilderBase(queue, @@ -484,10 +484,10 @@ class MatrixBlockBuilder(MatrixBlockBuilderBase): return mat -class P2PMatrixBlockBuilder(MatrixBlockBuilderBase): +class FarFieldBlockBuilder(MatrixBlockBuilderBase): def __init__(self, queue, dep_expr, other_dep_exprs, dep_source, places, context, index_set, exclude_self=True): - super(P2PMatrixBlockBuilder, self).__init__(queue, + super(FarFieldBlockBuilder, self).__init__(queue, dep_expr, other_dep_exprs, dep_source, places, context, index_set) self.dummy = MatrixBlockBuilderBase(queue, diff --git a/test/test_matrix.py b/test/test_matrix.py index 3becd048..aac9b86e 100644 --- a/test/test_matrix.py +++ b/test/test_matrix.py @@ -203,8 +203,8 @@ def test_p2p_block_builder(ctx_factory, factor, ndim, lpot_id, from sumpy.tools import MatrixBlockIndexRanges index_set = MatrixBlockIndexRanges(ctx, tgtindices, srcindices) - from pytential.symbolic.matrix import P2PMatrixBlockBuilder - mbuilder = P2PMatrixBlockBuilder(queue, + from pytential.symbolic.matrix import FarFieldBlockBuilder + mbuilder = FarFieldBlockBuilder(queue, dep_expr=u_sym, other_dep_exprs=[], dep_source=places[DEFAULT_SOURCE], @@ -237,7 +237,7 @@ def test_p2p_block_builder(ctx_factory, factor, ndim, lpot_id, _, (ax1, ax2) = mp.subplots(1, 2, figsize=(10, 8), dpi=300, constrained_layout=True) ax1.imshow(blk_full) - ax1.set_title('P2PMatrixBlockBuilder') + ax1.set_title('FarFieldBlockBuilder') ax2.imshow(mat_full) ax2.set_title('P2PMatrixBuilder') mp.savefig("test_p2p_block_{}d_{:.1f}.png".format(ndim, factor)) @@ -277,8 +277,8 @@ def test_qbx_block_builder(ctx_factory, ndim, lpot_id, visualize=False): from sumpy.tools import MatrixBlockIndexRanges index_set = MatrixBlockIndexRanges(ctx, tgtindices, srcindices) - from pytential.symbolic.matrix import MatrixBlockBuilder - mbuilder = MatrixBlockBuilder(queue, + from pytential.symbolic.matrix import NearFieldBlockBuilder + mbuilder = NearFieldBlockBuilder(queue, dep_expr=u_sym, other_dep_exprs=[], dep_source=places[DEFAULT_SOURCE], @@ -313,7 +313,7 @@ def test_qbx_block_builder(ctx_factory, ndim, lpot_id, visualize=False): ax1.imshow(mat_full) ax1.set_title('MatrixBuilder') ax2.imshow(blk_full) - ax2.set_title('MatrixBlockBuilder') + ax2.set_title('NearFieldBlockBuilder') mp.savefig("test_qbx_block_builder.png", dpi=300) for i in range(index_set.nblocks): -- GitLab From cf2785b775360fd2322c60112b7909d45e391fd5 Mon Sep 17 00:00:00 2001 From: Alex Fikl Date: Thu, 26 Jul 2018 17:54:48 -0500 Subject: [PATCH 09/34] matrix: allow building matrix for quad_stage2_density_discr --- pytential/symbolic/execution.py | 172 ++++++++++++++++++------------- pytential/symbolic/mappers.py | 4 + pytential/symbolic/matrix.py | 23 ++--- pytential/symbolic/primitives.py | 29 +++++- test/test_matrix.py | 46 ++++++++- 5 files changed, 187 insertions(+), 87 deletions(-) diff --git a/pytential/symbolic/execution.py b/pytential/symbolic/execution.py index f658a6aa..d5ab2c38 100644 --- a/pytential/symbolic/execution.py +++ b/pytential/symbolic/execution.py @@ -36,6 +36,7 @@ import pyopencl.clmath # noqa from loopy.version import MOST_RECENT_LANGUAGE_VERSION from pytools import memoize_in +from pytential.symbolic.primitives import DEFAULT_SOURCE, DEFAULT_TARGET # FIXME caches: fix up queues @@ -286,11 +287,11 @@ def _domains_default(nresults, places, domains, default_val): raise RuntimeError("'domains is None' requires " "default domain to be defined") dom_name = default_val - return nresults*[dom_name] + return nresults * [dom_name] elif not isinstance(domains, (list, tuple)): dom_name = domains - return nresults*[dom_name] + return nresults * [dom_name] else: return domains @@ -300,31 +301,43 @@ def _domains_default(nresults, places, domains, default_val): # {{{ bound expression +def _get_discretization(places, where, default_discr="density_discr"): + from pytential.source import LayerPotentialSourceBase + from pytential.symbolic.primitives import ( + _QBXSource, _QBXSourceStage2, _QBXSourceQuadStage2) + + if isinstance(where, _QBXSourceStage2): + name = "stage2_density_discr" + elif isinstance(where, _QBXSourceQuadStage2): + name = "quad_stage2_density_discr" + else: + name = default_discr + + if isinstance(where, _QBXSource): + where = where.where + discr = places[where] + + if isinstance(discr, LayerPotentialSourceBase): + return discr, getattr(discr, name) + return discr, discr + + class BoundExpression: - def __init__(self, optemplate, places): - self.optemplate = optemplate + def __init__(self, places, sym_op_expr, sym_op_args=None): self.places = places + self.sym_op_expr = sym_op_expr + self.sym_op_args = sym_op_args self.caches = {} from pytential.symbolic.compiler import OperatorCompiler - self.code = OperatorCompiler(self.places)(optemplate) + self.code = OperatorCompiler(self.places)(sym_op_expr) def get_cache(self, name): return self.caches.setdefault(name, {}) def get_discretization(self, where): - from pytential.symbolic.primitives import _QBXSourceStage2 - if isinstance(where, _QBXSourceStage2): - lpot_source = self.places[where.where] - return lpot_source.stage2_density_discr - - discr = self.places[where] - - from pytential.source import LayerPotentialSourceBase - if isinstance(discr, LayerPotentialSourceBase): - discr = discr.density_discr - + _, discr = _get_discretization(self.places, where) return discr def scipy_op(self, queue, arg_name, dtype, domains=None, **extra_args): @@ -343,7 +356,6 @@ class BoundExpression: else: nresults = 1 - from pytential.symbolic.primitives import DEFAULT_TARGET domains = _domains_default(nresults, self.places, domains, DEFAULT_TARGET) @@ -375,16 +387,32 @@ class BoundExpression: # {{{ expression prep -def prepare_places(places): - from pytential.symbolic.primitives import DEFAULT_SOURCE, DEFAULT_TARGET +def _where_default(places, auto_where=None): + if auto_where is None: + if isinstance(places, dict) and DEFAULT_TARGET not in places: + return (DEFAULT_SOURCE, DEFAULT_SOURCE) + else: + return (DEFAULT_SOURCE, DEFAULT_TARGET) + + return auto_where + + +def prepare_places(places, auto_where=None): from meshmode.discretization import Discretization from pytential.source import LayerPotentialSourceBase from pytential.target import TargetBase + where_source, where_target = _where_default(places, auto_where) if isinstance(places, LayerPotentialSourceBase): + from pytential.symbolic.primitives import _QBXSourceQuadStage2 + if isinstance(where_target, _QBXSourceQuadStage2): + target_discr = places.quad_stage2_density_discr + else: + target_discr = places.density_discr + places = { DEFAULT_SOURCE: places, - DEFAULT_TARGET: places.density_discr, + DEFAULT_TARGET: target_discr, } elif isinstance(places, (Discretization, TargetBase)): places = { @@ -420,13 +448,9 @@ def prepare_expr(places, expr, auto_where=None): evaluations, find 'where' attributes automatically. """ - from pytential.symbolic.primitives import DEFAULT_SOURCE, DEFAULT_TARGET from pytential.source import LayerPotentialSourceBase - from pytential.symbolic.mappers import ( - ToTargetTagger, - DerivativeBinder, - ) + ToTargetTagger, DerivativeBinder) if auto_where is None: if DEFAULT_TARGET in places: @@ -459,12 +483,33 @@ def bind(places, expr, auto_where=None): places = prepare_places(places) expr = prepare_expr(places, expr, auto_where=auto_where) - return BoundExpression(expr, places) + return BoundExpression(places, expr) # {{{ matrix building -def build_matrix(queue, places, expr, input_exprs, domains=None, +def prepare_expression(places, exprs, input_exprs, + domains=None, auto_where=None): + auto_where = _where_default(places, auto_where) + places = prepare_places(places, auto_where=auto_where) + exprs = prepare_expr(places, exprs, auto_where=auto_where) + + from pytools.obj_array import is_obj_array, make_obj_array + if not is_obj_array(exprs): + exprs = make_obj_array([exprs]) + try: + input_exprs = list(input_exprs) + except TypeError: + # not iterable, wrap in a list + input_exprs = [input_exprs] + + domains = _domains_default(len(input_exprs), places, domains, + DEFAULT_SOURCE) + + return places, exprs, input_exprs, domains + + +def build_matrix(queue, places, exprs, input_exprs, domains=None, auto_where=None, context=None): """ :arg queue: a :class:`pyopencl.CommandQueue` used to synchronize @@ -486,62 +531,49 @@ def build_matrix(queue, places, expr, input_exprs, domains=None, evaluations, find 'where' attributes automatically. """ + from pytools import single_valued + from pytential.symbolic.matrix import MatrixBuilder, is_zero + if context is None: context = {} - places = prepare_places(places) - expr = prepare_expr(places, expr, auto_where=auto_where) - - from pytools.obj_array import is_obj_array, make_obj_array - if not is_obj_array(expr): - expr = make_obj_array([expr]) - try: - input_exprs = list(input_exprs) - except TypeError: - # not iterable, wrap in a list - input_exprs = [input_exprs] - - from pytential.symbolic.primitives import DEFAULT_SOURCE - domains = _domains_default(len(input_exprs), places, domains, - DEFAULT_SOURCE) - - nblock_rows = len(expr) + auto_where = _where_default(places, auto_where) + places, exprs, input_exprs, domains = \ + prepare_expression(places, exprs, input_exprs, + domains=domains, + auto_where=auto_where) + auto_where = { + DEFAULT_SOURCE: auto_where[0], + DEFAULT_TARGET: auto_where[1] + } + + nblock_rows = len(exprs) nblock_columns = len(input_exprs) - blocks = np.zeros((nblock_rows, nblock_columns), dtype=np.object) - from pytential.symbolic.matrix import MatrixBuilder, is_zero - dtypes = [] - for ibcol in range(nblock_columns): + dep_source, dep_discr = \ + _get_discretization(places, auto_where[domains[ibcol]]) + mbuilder = MatrixBuilder( queue, dep_expr=input_exprs[ibcol], - other_dep_exprs=( - input_exprs[:ibcol] - + - input_exprs[ibcol+1:]), - dep_source=places[domains[ibcol]], + other_dep_exprs=(input_exprs[:ibcol] + + input_exprs[ibcol + 1:]), + dep_source=dep_source, + dep_discr=dep_discr, places=places, context=context) for ibrow in range(nblock_rows): - block = mbuilder(expr[ibrow]) - - assert ( - is_zero(block) - or isinstance(block, np.ndarray)) - if isinstance(block, np.ndarray): - dtypes.append(block.dtype) + block = mbuilder(exprs[ibrow]) + assert is_zero(block) or isinstance(block, np.ndarray) blocks[ibrow, ibcol] = block - - if isinstance(block, cl.array.Array): + if isinstance(block, np.ndarray): dtypes.append(block.dtype) - from pytools import single_valued - block_row_counts = [ single_valued( blocks[ibrow, ibcol].shape[0] @@ -550,20 +582,20 @@ def build_matrix(queue, places, expr, input_exprs, domains=None, for ibrow in range(nblock_rows)] block_col_counts = [ - places[domains[ibcol]].density_discr.nnodes + single_valued( + blocks[ibrow, ibcol].shape[1] + for ibrow in range(nblock_rows) + if not is_zero(blocks[ibrow, ibcol])) for ibcol in range(nblock_columns)] # "block row starts"/"block column starts" brs = np.cumsum([0] + block_row_counts) bcs = np.cumsum([0] + block_col_counts) - result = np.zeros( - (sum(block_row_counts), sum(block_col_counts)), - dtype=np.find_common_type(dtypes, [])) - + result = np.zeros((brs[-1], bcs[-1]), dtype=np.find_common_type(dtypes, [])) for ibcol in range(nblock_columns): for ibrow in range(nblock_rows): - result[brs[ibrow]:brs[ibrow+1], bcs[ibcol]:bcs[ibcol+1]] = \ + result[brs[ibrow]:brs[ibrow + 1], bcs[ibcol]:bcs[ibcol + 1]] = \ blocks[ibrow, ibcol] return cl.array.to_device(queue, result) diff --git a/pytential/symbolic/mappers.py b/pytential/symbolic/mappers.py index 9fdbaae4..69349479 100644 --- a/pytential/symbolic/mappers.py +++ b/pytential/symbolic/mappers.py @@ -453,8 +453,12 @@ class QBXPreprocessor(IdentityMapper): # {{{ stringifier def stringify_where(where): + if isinstance(where, prim._QBXSourceStage1): + return "stage1(%s)" % stringify_where(where.where) if isinstance(where, prim._QBXSourceStage2): return "stage2(%s)" % stringify_where(where.where) + if isinstance(where, prim._QBXSourceQuadStage2): + return "quad_stage2(%s)" % stringify_where(where.where) if where is None: return "?" diff --git a/pytential/symbolic/matrix.py b/pytential/symbolic/matrix.py index 4ea90d11..735782a8 100644 --- a/pytential/symbolic/matrix.py +++ b/pytential/symbolic/matrix.py @@ -92,15 +92,15 @@ def _get_kernel_args(mapper, kernel, expr, source): # We'll cheat and build the matrix on the host. class MatrixBuilder(EvaluationMapperBase): - def __init__(self, queue, dep_expr, other_dep_exprs, dep_source, places, - context): + def __init__(self, queue, dep_expr, other_dep_exprs, dep_source, dep_discr, + places, context): super(MatrixBuilder, self).__init__(context=context) self.queue = queue self.dep_expr = dep_expr self.other_dep_exprs = other_dep_exprs self.dep_source = dep_source - self.dep_discr = dep_source.density_discr + self.dep_discr = dep_discr self.places = places def map_variable(self, expr): @@ -188,11 +188,10 @@ class MatrixBuilder(EvaluationMapperBase): return vecs_and_scalars def map_int_g(self, expr): - source = self.places[expr.source] - target_discr = self.places[expr.target] - - if source.density_discr is not target_discr: - raise NotImplementedError() + from pytential.symbolic.execution import _get_discretization + source, source_discr = _get_discretization(self.places, expr.source, + default_discr="quad_stage2_density_discr") + _, target_discr = _get_discretization(self.places, expr.target) rec_density = self.rec(expr.density) if is_zero(rec_density): @@ -217,7 +216,7 @@ class MatrixBuilder(EvaluationMapperBase): _, (mat,) = mat_gen(self.queue, target_discr.nodes(), - source.quad_stage2_density_discr.nodes(), + source_discr.nodes(), get_centers_on_side(source, expr.qbx_forced_limit), expansion_radii=self.dep_source._expansion_radii("nsources"), **kernel_args) @@ -277,10 +276,10 @@ class MatrixBuilder(EvaluationMapperBase): # {{{ p2p matrix builder class P2PMatrixBuilder(MatrixBuilder): - def __init__(self, queue, dep_expr, other_dep_exprs, dep_source, places, - context, exclude_self=True): + def __init__(self, queue, dep_expr, other_dep_exprs, dep_source, dep_discr, + places, context, exclude_self=True): super(P2PMatrixBuilder, self).__init__(queue, - dep_expr, other_dep_exprs, dep_source, places, context) + dep_expr, other_dep_exprs, dep_source, dep_discr, places, context) self.exclude_self = exclude_self diff --git a/pytential/symbolic/primitives.py b/pytential/symbolic/primitives.py index 3a297fd1..909dbe5a 100644 --- a/pytential/symbolic/primitives.py +++ b/pytential/symbolic/primitives.py @@ -201,10 +201,10 @@ class DEFAULT_TARGET: # noqa pass -class _QBXSourceStage2(object): - """A symbolic 'where' specifier for the - :attr:`pytential.qbx.QBXLayerPotentialSource.stage2_density_discr` - of the layer potential source identified by :attr:`where`. +class _QBXSource(object): + """A symbolic 'where' specifier for the a density of a + :attr:`pytential.qbx.QBXLayerPotentialSource` + layer potential source identified by :attr:`where`. .. attribute:: where @@ -229,6 +229,27 @@ class _QBXSourceStage2(object): def __ne__(self, other): return not self.__eq__(other) + +class _QBXSourceStage1(_QBXSource): + """An explicit symbolic 'where' specifier for the + :attr:`pytential.qbx.QBXLayerPotentialSource.density_discr` + of the layer potential source identified by :attr:`where`. + """ + + +class _QBXSourceStage2(_QBXSource): + """A symbolic 'where' specifier for the + :attr:`pytential.qbx.QBXLayerPotentialSource.stage2_density_discr` + of the layer potential source identified by :attr:`where`. + """ + + +class _QBXSourceQuadStage2(_QBXSource): + """A symbolic 'where' specifier for the + :attr:`pytential.qbx.QBXLayerPotentialSource.quad_stage2_density_discr` + of the layer potential source identified by :attr:`where`. + """ + # }}} diff --git a/test/test_matrix.py b/test/test_matrix.py index aac9b86e..e168cbc3 100644 --- a/test/test_matrix.py +++ b/test/test_matrix.py @@ -38,7 +38,8 @@ from meshmode.mesh.generation import \ ellipse, NArmedStarfish, make_curve_mesh from pytential import bind, sym -from pytential.symbolic.primitives import DEFAULT_SOURCE +from pytential.symbolic.primitives import DEFAULT_SOURCE, DEFAULT_TARGET +from pytential.symbolic.primitives import _QBXSourceStage1, _QBXSourceQuadStage2 import pytest from pyopencl.tools import ( # noqa @@ -218,6 +219,7 @@ def test_p2p_block_builder(ctx_factory, factor, ndim, lpot_id, dep_expr=u_sym, other_dep_exprs=[], dep_source=places[DEFAULT_SOURCE], + dep_discr=places[DEFAULT_SOURCE].density_discr, places=places, context={}) mat = mbuilder(expr) @@ -292,6 +294,7 @@ def test_qbx_block_builder(ctx_factory, ndim, lpot_id, visualize=False): dep_expr=u_sym, other_dep_exprs=[], dep_source=places[DEFAULT_SOURCE], + dep_discr=places[DEFAULT_SOURCE].density_discr, places=places, context={}) mat = mbuilder(expr) @@ -326,6 +329,47 @@ def test_qbx_block_builder(ctx_factory, ndim, lpot_id, visualize=False): assert error < eps +@pytest.mark.parametrize('where', + [None, + (DEFAULT_SOURCE, DEFAULT_TARGET), + (_QBXSourceStage1(DEFAULT_SOURCE), + _QBXSourceStage1(DEFAULT_TARGET)), + (_QBXSourceQuadStage2(DEFAULT_SOURCE), + _QBXSourceQuadStage2(DEFAULT_TARGET))]) +def test_build_matrix_where(ctx_factory, where): + pytest.skip("wip") + ctx = ctx_factory() + queue = cl.CommandQueue(ctx) + + # prevent cache explosion + from sympy.core.cache import clear_cache + clear_cache() + + from test_linalg_proxy import _build_qbx_discr + qbx = _build_qbx_discr(queue, target_order=7, ndim=2) + op, u_sym, _ = _build_op(lpot_id=1, ndim=2) + + from pytential.symbolic.execution import build_matrix + mat = build_matrix(queue, qbx, op, u_sym, auto_where=where) + + if where is None: + source_where, target_where = DEFAULT_SOURCE, DEFAULT_TARGET + else: + source_where, target_where = where + + if isinstance(source_where, _QBXSourceQuadStage2): + n = qbx.quad_stage2_density_discr.nnodes + else: + n = qbx.density_discr.nnodes + + if isinstance(target_where, _QBXSourceQuadStage2): + m = qbx.quad_stage2_density_discr.nnodes + else: + m = qbx.density_discr.nnodes + + assert mat.shape == (m, n) + + if __name__ == "__main__": import sys if len(sys.argv) > 1: -- GitLab From bf8a02a38ff86261ffbf43fb9ac35d06a54f64b6 Mon Sep 17 00:00:00 2001 From: Alex Fikl Date: Sun, 29 Jul 2018 20:43:56 -0500 Subject: [PATCH 10/34] matrix: more work towards supporting different source and target where locations --- pytential/qbx/__init__.py | 2 +- pytential/symbolic/execution.py | 254 +++++++++++++++---------------- pytential/symbolic/mappers.py | 12 +- pytential/symbolic/matrix.py | 43 +++++- pytential/symbolic/primitives.py | 6 +- test/test_matrix.py | 15 +- 6 files changed, 176 insertions(+), 156 deletions(-) diff --git a/pytential/qbx/__init__.py b/pytential/qbx/__init__.py index 869f9222..2bbe3d36 100644 --- a/pytential/qbx/__init__.py +++ b/pytential/qbx/__init__.py @@ -564,7 +564,7 @@ class QBXLayerPotentialSource(LayerPotentialSourceBase): maxstretch = bind( self, sym._simplex_mapping_max_stretch_factor( self.ambient_dim, - where=sym._QBXSourceStage2(sym.DEFAULT_SOURCE)) + where=sym.QBXSourceStage2(sym.DEFAULT_SOURCE)) )(queue) maxstretch = utils.to_last_dim_length( self.stage2_density_discr, maxstretch, last_dim_length) diff --git a/pytential/symbolic/execution.py b/pytential/symbolic/execution.py index d5ab2c38..2ec5b05f 100644 --- a/pytential/symbolic/execution.py +++ b/pytential/symbolic/execution.py @@ -37,6 +37,8 @@ from loopy.version import MOST_RECENT_LANGUAGE_VERSION from pytools import memoize_in from pytential.symbolic.primitives import DEFAULT_SOURCE, DEFAULT_TARGET +from pytential.symbolic.primitives import ( + QBXSourceStage1, QBXSourceStage2, QBXSourceQuadStage2) # FIXME caches: fix up queues @@ -279,7 +281,7 @@ class MatVecOp: # }}} -# {{{ default for 'domains' parameter +# {{{ expression prep def _domains_default(nresults, places, domains, default_val): if domains is None: @@ -296,30 +298,130 @@ def _domains_default(nresults, places, domains, default_val): else: return domains + +def _where_default(places, auto_where=None): + if auto_where is None: + if not isinstance(places, dict): + return DEFAULT_SOURCE, DEFAULT_TARGET + if DEFAULT_TARGET in places: + return DEFAULT_SOURCE, DEFAULT_TARGET + return tuple(places.keys()) + + return auto_where + + +def prepare_places(places, auto_where=None): + from meshmode.discretization import Discretization + from pytential.source import LayerPotentialSourceBase + from pytential.target import TargetBase + + where_source, where_target = _where_default(places, auto_where=auto_where) + if isinstance(places, LayerPotentialSourceBase): + _, target_discr = _get_discretization(places, where_target) + places = { + where_source: places, + where_target: target_discr, + } + elif isinstance(places, (Discretization, TargetBase)): + places = { + where_target: places, + } + + elif isinstance(places, tuple): + source_discr, target_discr = places + places = { + where_source: source_discr, + where_target: target_discr, + } + del source_discr + del target_discr + + def cast_to_place(discr): + from pytential.target import TargetBase + from pytential.source import PotentialSource + if not isinstance(discr, (Discretization, TargetBase, PotentialSource)): + raise TypeError("must pass discretizations, " + "layer potential sources or targets as 'places'") + return discr + + return dict( + (key, cast_to_place(value)) + for key, value in six.iteritems(places)) + + +def prepare_expr(places, expr, auto_where=None): + """ + :arg places: result of :func:`prepare_places` + :arg auto_where: For simple source-to-self or source-to-target + evaluations, find 'where' attributes automatically. + """ + + from pytential.source import LayerPotentialSourceBase + from pytential.symbolic.mappers import ( + ToTargetTagger, DerivativeBinder) + + auto_where = _where_default(places, auto_where=auto_where) + if auto_where: + expr = ToTargetTagger(*auto_where)(expr) + + expr = DerivativeBinder()(expr) + + for name, place in six.iteritems(places): + if isinstance(place, LayerPotentialSourceBase): + expr = place.preprocess_optemplate(name, places, expr) + + return expr + + +def prepare_expression(places, exprs, input_exprs, + domains=None, auto_where=None): + auto_where = _where_default(places, auto_where) + places = prepare_places(places, auto_where=auto_where) + exprs = prepare_expr(places, exprs, auto_where=auto_where) + + from pytools.obj_array import is_obj_array, make_obj_array + if not is_obj_array(exprs): + exprs = make_obj_array([exprs]) + try: + input_exprs = list(input_exprs) + except TypeError: + # not iterable, wrap in a list + input_exprs = [input_exprs] + + domains = _domains_default(len(input_exprs), places, domains, auto_where[0]) + + return places, exprs, input_exprs, domains + # }}} # {{{ bound expression -def _get_discretization(places, where, default_discr="density_discr"): +def _get_discretization(places, where, default_source=QBXSourceStage1): from pytential.source import LayerPotentialSourceBase - from pytential.symbolic.primitives import ( - _QBXSource, _QBXSourceStage2, _QBXSourceQuadStage2) - if isinstance(where, _QBXSourceStage2): - name = "stage2_density_discr" - elif isinstance(where, _QBXSourceQuadStage2): - name = "quad_stage2_density_discr" - else: - name = default_discr + if where is DEFAULT_SOURCE or where is DEFAULT_TARGET: + where = default_source(where) - if isinstance(where, _QBXSource): - where = where.where - discr = places[where] + if isinstance(places, LayerPotentialSourceBase): + lpot = places + else: + try: + lpot = places[where.where] + except KeyError: + lpot = places[where] + is_lpot = isinstance(lpot, LayerPotentialSourceBase) + + if isinstance(where, QBXSourceStage1): + discr = lpot.density_discr if is_lpot else lpot + elif isinstance(where, QBXSourceStage2): + discr = lpot.stage2_density_discr if is_lpot else lpot + elif isinstance(where, QBXSourceQuadStage2): + discr = lpot.quad_stage2_density_discr if is_lpot else lpot + else: + raise ValueError("Unknown 'where': {}".format(type(where))) - if isinstance(discr, LayerPotentialSourceBase): - return discr, getattr(discr, name) - return discr, discr + return lpot, discr class BoundExpression: @@ -382,95 +484,6 @@ class BoundExpression: exec_mapper = EvaluationMapper(self, queue, args) return self.code.execute(exec_mapper) -# }}} - - -# {{{ expression prep - -def _where_default(places, auto_where=None): - if auto_where is None: - if isinstance(places, dict) and DEFAULT_TARGET not in places: - return (DEFAULT_SOURCE, DEFAULT_SOURCE) - else: - return (DEFAULT_SOURCE, DEFAULT_TARGET) - - return auto_where - - -def prepare_places(places, auto_where=None): - from meshmode.discretization import Discretization - from pytential.source import LayerPotentialSourceBase - from pytential.target import TargetBase - - where_source, where_target = _where_default(places, auto_where) - if isinstance(places, LayerPotentialSourceBase): - from pytential.symbolic.primitives import _QBXSourceQuadStage2 - if isinstance(where_target, _QBXSourceQuadStage2): - target_discr = places.quad_stage2_density_discr - else: - target_discr = places.density_discr - - places = { - DEFAULT_SOURCE: places, - DEFAULT_TARGET: target_discr, - } - elif isinstance(places, (Discretization, TargetBase)): - places = { - DEFAULT_TARGET: places, - } - - elif isinstance(places, tuple): - source_discr, target_discr = places - places = { - DEFAULT_SOURCE: source_discr, - DEFAULT_TARGET: target_discr, - } - del source_discr - del target_discr - - def cast_to_place(discr): - from pytential.target import TargetBase - from pytential.source import PotentialSource - if not isinstance(discr, (Discretization, TargetBase, PotentialSource)): - raise TypeError("must pass discretizations, " - "layer potential sources or targets as 'places'") - return discr - - return dict( - (key, cast_to_place(value)) - for key, value in six.iteritems(places)) - - -def prepare_expr(places, expr, auto_where=None): - """ - :arg places: result of :func:`prepare_places` - :arg auto_where: For simple source-to-self or source-to-target - evaluations, find 'where' attributes automatically. - """ - - from pytential.source import LayerPotentialSourceBase - from pytential.symbolic.mappers import ( - ToTargetTagger, DerivativeBinder) - - if auto_where is None: - if DEFAULT_TARGET in places: - auto_where = DEFAULT_SOURCE, DEFAULT_TARGET - else: - auto_where = DEFAULT_SOURCE, DEFAULT_SOURCE - - if auto_where: - expr = ToTargetTagger(*auto_where)(expr) - - expr = DerivativeBinder()(expr) - - for name, place in six.iteritems(places): - if isinstance(place, LayerPotentialSourceBase): - expr = place.preprocess_optemplate(name, places, expr) - - return expr - -# }}} - def bind(places, expr, auto_where=None): """ @@ -485,29 +498,10 @@ def bind(places, expr, auto_where=None): expr = prepare_expr(places, expr, auto_where=auto_where) return BoundExpression(places, expr) +# }}} -# {{{ matrix building - -def prepare_expression(places, exprs, input_exprs, - domains=None, auto_where=None): - auto_where = _where_default(places, auto_where) - places = prepare_places(places, auto_where=auto_where) - exprs = prepare_expr(places, exprs, auto_where=auto_where) - - from pytools.obj_array import is_obj_array, make_obj_array - if not is_obj_array(exprs): - exprs = make_obj_array([exprs]) - try: - input_exprs = list(input_exprs) - except TypeError: - # not iterable, wrap in a list - input_exprs = [input_exprs] - - domains = _domains_default(len(input_exprs), places, domains, - DEFAULT_SOURCE) - - return places, exprs, input_exprs, domains +# {{{ matrix building def build_matrix(queue, places, exprs, input_exprs, domains=None, auto_where=None, context=None): @@ -537,15 +531,11 @@ def build_matrix(queue, places, exprs, input_exprs, domains=None, if context is None: context = {} - auto_where = _where_default(places, auto_where) + auto_where = _where_default(places, auto_where=auto_where) places, exprs, input_exprs, domains = \ prepare_expression(places, exprs, input_exprs, domains=domains, auto_where=auto_where) - auto_where = { - DEFAULT_SOURCE: auto_where[0], - DEFAULT_TARGET: auto_where[1] - } nblock_rows = len(exprs) nblock_columns = len(input_exprs) @@ -554,7 +544,7 @@ def build_matrix(queue, places, exprs, input_exprs, domains=None, dtypes = [] for ibcol in range(nblock_columns): dep_source, dep_discr = \ - _get_discretization(places, auto_where[domains[ibcol]]) + _get_discretization(places, domains[ibcol]) mbuilder = MatrixBuilder( queue, diff --git a/pytential/symbolic/mappers.py b/pytential/symbolic/mappers.py index 69349479..45f3b84f 100644 --- a/pytential/symbolic/mappers.py +++ b/pytential/symbolic/mappers.py @@ -292,8 +292,10 @@ class ToTargetTagger(LocationTagger): """ def __init__(self, default_source, default_target): - LocationTagger.__init__(self, default_target) - self.operand_rec = LocationTagger(default_source) + LocationTagger.__init__(self, default_target, + default_source=default_source) + self.operand_rec = LocationTagger(default_source, + default_source=default_source) # }}} @@ -453,11 +455,11 @@ class QBXPreprocessor(IdentityMapper): # {{{ stringifier def stringify_where(where): - if isinstance(where, prim._QBXSourceStage1): + if isinstance(where, prim.QBXSourceStage1): return "stage1(%s)" % stringify_where(where.where) - if isinstance(where, prim._QBXSourceStage2): + if isinstance(where, prim.QBXSourceStage2): return "stage2(%s)" % stringify_where(where.where) - if isinstance(where, prim._QBXSourceQuadStage2): + if isinstance(where, prim.QBXSourceQuadStage2): return "quad_stage2(%s)" % stringify_where(where.where) if where is None: diff --git a/pytential/symbolic/matrix.py b/pytential/symbolic/matrix.py index 735782a8..d6f15922 100644 --- a/pytential/symbolic/matrix.py +++ b/pytential/symbolic/matrix.py @@ -83,6 +83,30 @@ def _get_kernel_args(mapper, kernel, expr, source): return kernel_args + +def _weights_and_area_elements(queue, source, where): + if isinstance(where, sym.QBXSourceQuadStage2): + return source.weights_and_area_elements() + + # NOTE: copied from `weights_and_area_elements`, but using the + # discretization given by `where` and no interpolation + from pytential.symbolic.execution import _get_discretization + discr = _get_discretization(source, where) + + area = bind(discr, sym.area_element(source.ambient_dim, source.dim))(queue) + qweight = bind(discr, sym.QWeight())(queue) + + return area * qweight + + +def _get_centers_on_side(queue, source, qbx_forced_limit): + pass + + +def _get_expansion_radii(queue, source): + pass + + # }}} @@ -189,8 +213,12 @@ class MatrixBuilder(EvaluationMapperBase): def map_int_g(self, expr): from pytential.symbolic.execution import _get_discretization - source, source_discr = _get_discretization(self.places, expr.source, - default_discr="quad_stage2_density_discr") + if expr.source is sym.DEFAULT_SOURCE: + where_source = sym.QBXSourceQuadStage2(expr.source) + else: + where_source = expr.source + + source, source_discr = _get_discretization(self.places, where_source) _, target_discr = _get_discretization(self.places, expr.target) rec_density = self.rec(expr.density) @@ -222,13 +250,14 @@ class MatrixBuilder(EvaluationMapperBase): **kernel_args) mat = mat.get() - waa = source.weights_and_area_elements().get(queue=self.queue) - mat[:, :] *= waa + waa = _weights_and_area_elements(self.queue, source, where_source) + mat[:, :] *= waa.get(self.queue) - resampler = source.direct_resampler - resample_mat = resampler.full_resample_matrix(self.queue).get(self.queue) + if target_discr.nnodes != source_discr.nnodes: + resampler = source.direct_resampler + resample_mat = resampler.full_resample_matrix(self.queue).get(self.queue) + mat = mat.dot(resample_mat) - mat = mat.dot(resample_mat) mat = mat.dot(rec_density) return mat diff --git a/pytential/symbolic/primitives.py b/pytential/symbolic/primitives.py index 909dbe5a..2288ab85 100644 --- a/pytential/symbolic/primitives.py +++ b/pytential/symbolic/primitives.py @@ -230,21 +230,21 @@ class _QBXSource(object): return not self.__eq__(other) -class _QBXSourceStage1(_QBXSource): +class QBXSourceStage1(_QBXSource): """An explicit symbolic 'where' specifier for the :attr:`pytential.qbx.QBXLayerPotentialSource.density_discr` of the layer potential source identified by :attr:`where`. """ -class _QBXSourceStage2(_QBXSource): +class QBXSourceStage2(_QBXSource): """A symbolic 'where' specifier for the :attr:`pytential.qbx.QBXLayerPotentialSource.stage2_density_discr` of the layer potential source identified by :attr:`where`. """ -class _QBXSourceQuadStage2(_QBXSource): +class QBXSourceQuadStage2(_QBXSource): """A symbolic 'where' specifier for the :attr:`pytential.qbx.QBXLayerPotentialSource.quad_stage2_density_discr` of the layer potential source identified by :attr:`where`. diff --git a/test/test_matrix.py b/test/test_matrix.py index e168cbc3..e3cc0c8d 100644 --- a/test/test_matrix.py +++ b/test/test_matrix.py @@ -39,7 +39,7 @@ from meshmode.mesh.generation import \ from pytential import bind, sym from pytential.symbolic.primitives import DEFAULT_SOURCE, DEFAULT_TARGET -from pytential.symbolic.primitives import _QBXSourceStage1, _QBXSourceQuadStage2 +from pytential.symbolic.primitives import QBXSourceStage1, QBXSourceQuadStage2 import pytest from pyopencl.tools import ( # noqa @@ -332,12 +332,11 @@ def test_qbx_block_builder(ctx_factory, ndim, lpot_id, visualize=False): @pytest.mark.parametrize('where', [None, (DEFAULT_SOURCE, DEFAULT_TARGET), - (_QBXSourceStage1(DEFAULT_SOURCE), - _QBXSourceStage1(DEFAULT_TARGET)), - (_QBXSourceQuadStage2(DEFAULT_SOURCE), - _QBXSourceQuadStage2(DEFAULT_TARGET))]) + (QBXSourceStage1(DEFAULT_SOURCE), + QBXSourceStage1(DEFAULT_TARGET)), + (QBXSourceQuadStage2(DEFAULT_SOURCE), + QBXSourceQuadStage2(DEFAULT_TARGET))]) def test_build_matrix_where(ctx_factory, where): - pytest.skip("wip") ctx = ctx_factory() queue = cl.CommandQueue(ctx) @@ -357,12 +356,12 @@ def test_build_matrix_where(ctx_factory, where): else: source_where, target_where = where - if isinstance(source_where, _QBXSourceQuadStage2): + if isinstance(source_where, QBXSourceQuadStage2): n = qbx.quad_stage2_density_discr.nnodes else: n = qbx.density_discr.nnodes - if isinstance(target_where, _QBXSourceQuadStage2): + if isinstance(target_where, QBXSourceQuadStage2): m = qbx.quad_stage2_density_discr.nnodes else: m = qbx.density_discr.nnodes -- GitLab From 72e0821d84e074692bf988ae45ed6e3c8e7c9e59 Mon Sep 17 00:00:00 2001 From: Alex Fikl Date: Tue, 31 Jul 2018 15:07:22 -0500 Subject: [PATCH 11/34] matrix: add code to associate stage2 targets to expansion centers --- pytential/symbolic/mappers.py | 11 +++----- pytential/symbolic/matrix.py | 47 ++++++++++++++++++++++++----------- test/test_matrix.py | 44 ++++++++++++++++++++++++++++++-- 3 files changed, 79 insertions(+), 23 deletions(-) diff --git a/pytential/symbolic/mappers.py b/pytential/symbolic/mappers.py index 45f3b84f..12b87894 100644 --- a/pytential/symbolic/mappers.py +++ b/pytential/symbolic/mappers.py @@ -399,19 +399,16 @@ class QBXPreprocessor(IdentityMapper): self.places = places def map_int_g(self, expr): - source = self.places[self.source_name] - target_discr = self.places[expr.target] - - from pytential.source import LayerPotentialSourceBase - if isinstance(target_discr, LayerPotentialSourceBase): - target_discr = target_discr.density_discr + from pytential.symbolic.execution import _get_discretization + source, source_discr = _get_discretization(self.places, self.source_name) + _, target_discr = _get_discretization(self.places, expr.target) if expr.qbx_forced_limit == 0: raise ValueError("qbx_forced_limit == 0 was a bad idea and " "is no longer supported. Use qbx_forced_limit == 'avg' " "to request two-sided averaging explicitly if needed.") - is_self = source.density_discr is target_discr + is_self = source_discr is target_discr expr = expr.copy( kernel=expr.kernel, diff --git a/pytential/symbolic/matrix.py b/pytential/symbolic/matrix.py index d6f15922..356dacc6 100644 --- a/pytential/symbolic/matrix.py +++ b/pytential/symbolic/matrix.py @@ -99,13 +99,26 @@ def _weights_and_area_elements(queue, source, where): return area * qweight -def _get_centers_on_side(queue, source, qbx_forced_limit): - pass - - -def _get_expansion_radii(queue, source): - pass - +def _get_centers_and_expansion_radii(queue, source, target_discr, qbx_forced_limit): + from pytential.qbx.utils import get_interleaved_centers + centers = get_interleaved_centers(queue, source) + radii = source._expansion_radii('nsources') + + from pytential.qbx.target_assoc import associate_targets_to_qbx_centers + code_container = source.target_association_code_container + assoc = associate_targets_to_qbx_centers( + source, + code_container.get_wrangler(queue), + [(target_discr, qbx_forced_limit)], + target_association_tolerance=1.0e-1) + + centers = [cl.array.take(c, assoc.target_to_center, queue=queue) + for c in centers] + radii = cl.array.take(radii, + (assoc.target_to_center.with_queue(queue) / 2.0).astype(np.int), + queue=queue) + + return centers, radii # }}} @@ -212,14 +225,19 @@ class MatrixBuilder(EvaluationMapperBase): return vecs_and_scalars def map_int_g(self, expr): - from pytential.symbolic.execution import _get_discretization if expr.source is sym.DEFAULT_SOURCE: where_source = sym.QBXSourceQuadStage2(expr.source) else: where_source = expr.source + if expr.target is sym.DEFAULT_TARGET: + where_target = sym.QBXSourceStage1(expr.target) + else: + where_target = expr.target + + from pytential.symbolic.execution import _get_discretization source, source_discr = _get_discretization(self.places, where_source) - _, target_discr = _get_discretization(self.places, expr.target) + _, target_discr = _get_discretization(self.places, where_target) rec_density = self.rec(expr.density) if is_zero(rec_density): @@ -239,14 +257,15 @@ class MatrixBuilder(EvaluationMapperBase): mat_gen = LayerPotentialMatrixGenerator( self.queue.context, (local_expn,)) - from pytential.qbx.utils import get_centers_on_side assert abs(expr.qbx_forced_limit) > 0 + centers, radii = _get_centers_and_expansion_radii(self.queue, + source, target_discr, expr.qbx_forced_limit) _, (mat,) = mat_gen(self.queue, - target_discr.nodes(), - source_discr.nodes(), - get_centers_on_side(source, expr.qbx_forced_limit), - expansion_radii=self.dep_source._expansion_radii("nsources"), + targets=target_discr.nodes(), + sources=source_discr.nodes(), + centers=centers, + expansion_radii=radii, **kernel_args) mat = mat.get() diff --git a/test/test_matrix.py b/test/test_matrix.py index e3cc0c8d..b266c2e8 100644 --- a/test/test_matrix.py +++ b/test/test_matrix.py @@ -336,7 +336,7 @@ def test_qbx_block_builder(ctx_factory, ndim, lpot_id, visualize=False): QBXSourceStage1(DEFAULT_TARGET)), (QBXSourceQuadStage2(DEFAULT_SOURCE), QBXSourceQuadStage2(DEFAULT_TARGET))]) -def test_build_matrix_where(ctx_factory, where): +def test_build_matrix_where(ctx_factory, where, visualize=False): ctx = ctx_factory() queue = cl.CommandQueue(ctx) @@ -345,8 +345,48 @@ def test_build_matrix_where(ctx_factory, where): clear_cache() from test_linalg_proxy import _build_qbx_discr - qbx = _build_qbx_discr(queue, target_order=7, ndim=2) + qbx = _build_qbx_discr(queue, nelements=8, target_order=2, ndim=2, + curve_f=partial(ellipse, 1.0)) op, u_sym, _ = _build_op(lpot_id=1, ndim=2) + qbx_forced_limit = -1 + + # associate quad_stage2 targets to centers + from pytential.qbx.target_assoc import associate_targets_to_qbx_centers + + code_container = qbx.target_association_code_container + target_discr = qbx.quad_stage2_density_discr + target_assoc = associate_targets_to_qbx_centers( + qbx, + code_container.get_wrangler(queue), + [(target_discr, qbx_forced_limit)], + target_association_tolerance=1.0e-1).get(queue) + target_assoc = target_assoc.target_to_center + + if qbx.ambient_dim == 2 and visualize: + import matplotlib.pyplot as pt + from pytential.qbx.utils import get_interleaved_centers + sources = qbx.density_discr.nodes().get(queue) + targets = target_discr.nodes().get(queue) + centers = get_interleaved_centers(queue, qbx) + centers = np.vstack(c.get(queue) for c in centers) + radii = qbx._expansion_radii('nsources').get(queue) + + for i in range(centers[0].size): + itgt = np.where(target_assoc == i)[0] + if not len(itgt): + continue + + pt.figure(figsize=(10, 8), dpi=300) + pt.plot(sources[0], sources[1], 'k'); + pt.plot(targets[0], targets[1], 'ko'); + + line = pt.plot(targets[0, itgt], targets[1, itgt], 'o') + c = pt.Circle([centers[0][i], centers[1][i]], radii[i // 2], + color=line[0].get_color(), alpha=0.5) + pt.gca().add_artist(c) + + pt.savefig('test_assoc_quad_targets_to_centers_{:05}.png'.format(i // 2)) + pt.close() from pytential.symbolic.execution import build_matrix mat = build_matrix(queue, qbx, op, u_sym, auto_where=where) -- GitLab From bcebd37ec86531ba1cd58ffb7bf95a44d8040597 Mon Sep 17 00:00:00 2001 From: Alex Fikl Date: Tue, 31 Jul 2018 15:35:29 -0500 Subject: [PATCH 12/34] matrix: port other matrix builders to different where types --- pytential/symbolic/matrix.py | 126 ++++++++++++++++------------------- test/test_matrix.py | 23 ++----- 2 files changed, 64 insertions(+), 85 deletions(-) diff --git a/pytential/symbolic/matrix.py b/pytential/symbolic/matrix.py index 356dacc6..4e6ded89 100644 --- a/pytential/symbolic/matrix.py +++ b/pytential/symbolic/matrix.py @@ -86,7 +86,7 @@ def _get_kernel_args(mapper, kernel, expr, source): def _weights_and_area_elements(queue, source, where): if isinstance(where, sym.QBXSourceQuadStage2): - return source.weights_and_area_elements() + return source.weights_and_area_elements().with_queue(queue) # NOTE: copied from `weights_and_area_elements`, but using the # discretization given by `where` and no interpolation @@ -225,19 +225,18 @@ class MatrixBuilder(EvaluationMapperBase): return vecs_and_scalars def map_int_g(self, expr): - if expr.source is sym.DEFAULT_SOURCE: - where_source = sym.QBXSourceQuadStage2(expr.source) - else: - where_source = expr.source + where_source = expr.source + if where_source is DEFAULT_SOURCE: + where_source = sym.QBXSourceQuadStage2(where_source) - if expr.target is sym.DEFAULT_TARGET: + where_target = expr.target + if where_target is DEFAULT_TARGET where_target = sym.QBXSourceStage1(expr.target) - else: - where_target = expr.target from pytential.symbolic.execution import _get_discretization source, source_discr = _get_discretization(self.places, where_source) _, target_discr = _get_discretization(self.places, where_target) + assert source_discr.nnodes >= target_discr.nnodes rec_density = self.rec(expr.density) if is_zero(rec_density): @@ -332,11 +331,18 @@ class P2PMatrixBuilder(MatrixBuilder): self.exclude_self = exclude_self def map_int_g(self, expr): - source = self.places[expr.source] - target_discr = self.places[expr.target] + where_source = expr.source + if where_source is DEFAULT_SOURCE: + where_source = sym.QBXSourceQuadStage2(where_source) - if source.density_discr is not target_discr: - raise NotImplementedError() + where_target = expr.target + if where_target is DEFAULT_TARGET + where_target = sym.QBXSourceStage1(expr.target) + + from pytential.symbolic.execution import _get_discretization + source, source_discr = _get_discretization(self.places, where_source) + _, target_discr = _get_discretization(self.places, where_target) + assert source_discr.nnodes >= target_discr.nnodes rec_density = self.rec(expr.density) if is_zero(rec_density): @@ -346,11 +352,8 @@ class P2PMatrixBuilder(MatrixBuilder): if len(rec_density.shape) != 2: raise NotImplementedError("layer potentials on non-variables") - try: - kernel = expr.kernel.inner_kernel - except AttributeError: - kernel = expr.kernel - kernel_args = _get_kernel_args(self, kernel, expr, source) + kernel_args = _get_kernel_args(self, + kernel.get_base_kernel(), expr, source) if self.exclude_self: kernel_args["target_to_source"] = \ cl.array.arange(self.queue, 0, target_discr.nnodes, dtype=np.int) @@ -361,7 +364,7 @@ class P2PMatrixBuilder(MatrixBuilder): _, (mat,) = mat_gen(self.queue, targets=target_discr.nodes(), - sources=source.density_discr.nodes(), + sources=source_discr.nodes(), **kernel_args) mat = mat.get() @@ -451,17 +454,6 @@ class NearFieldBlockBuilder(MatrixBlockBuilderBase): self.dummy = MatrixBlockBuilderBase(queue, dep_expr, other_dep_exprs, dep_source, places, context, index_set) - @memoize_method - def _oversampled_index_set(self, resampler): - from pytential.linalg.proxy import partition_from_coarse - srcindices = partition_from_coarse(resampler, self.index_set) - - from sumpy.tools import MatrixBlockIndexRanges - ovsm_index_set = MatrixBlockIndexRanges(self.queue.context, - self.index_set.row, srcindices) - - return ovsm_index_set - def _map_dep_variable(self): tgtindices = self.index_set.row.indices.get(self.queue).reshape(-1, 1) srcindices = self.index_set.col.indices.get(self.queue).reshape(1, -1) @@ -469,10 +461,19 @@ class NearFieldBlockBuilder(MatrixBlockBuilderBase): return np.equal(tgtindices, srcindices).astype(np.float64) def map_int_g(self, expr): - source = self.places[expr.source] - target_discr = self.places[expr.target] + where_source = expr.source + if where_source is DEFAULT_SOURCE: + where_source = sym.QBXSourceQuadStage2(where_source) + + where_target = expr.target + if where_target is DEFAULT_TARGET + where_target = sym.QBXSourceQuadStage2(expr.target) + + from pytential.symbolic.execution import _get_discretization + source, source_discr = _get_discretization(self.places, where_source) + _, target_discr = _get_discretization(self.places, where_target) - if source.density_discr is not target_discr: + if source_discr is not target_discr: raise NotImplementedError() rec_density = self.dummy.rec(expr.density) @@ -483,9 +484,6 @@ class NearFieldBlockBuilder(MatrixBlockBuilderBase): if len(rec_density.shape) != 2: raise NotImplementedError("layer potentials on non-variables") - resampler = source.direct_resampler - ovsm_index_set = self._oversampled_index_set(resampler) - kernel = expr.kernel kernel_args = _get_layer_potential_args(self, expr, source) @@ -498,35 +496,21 @@ class NearFieldBlockBuilder(MatrixBlockBuilderBase): from pytential.qbx.utils import get_centers_on_side assert abs(expr.qbx_forced_limit) > 0 + centers, radii = _get_centers_and_expansion_radii(queue, + source, target_discr, expr.qbx_forced_limit) _, (mat,) = mat_gen(self.queue, - target_discr.nodes(), - source.quad_stage2_density_discr.nodes(), - get_centers_on_side(source, expr.qbx_forced_limit), - expansion_radii=self.dep_source._expansion_radii("nsources"), - index_set=ovsm_index_set, + targets=target_discr.nodes(), + sources=source_discr.nodes(), + centers=centers, + expansion_radii=radii, + index_set=index_set, **kernel_args) - waa = source.weights_and_area_elements().with_queue(self.queue) - mat *= waa[ovsm_index_set.linear_col_indices] - - from sumpy.tools import MatrixBlockIndexRanges - ovsm_col_index = ovsm_index_set.get(self.queue) - ovsm_row_index = MatrixBlockIndexRanges(ovsm_col_index.cl_context, - ovsm_col_index.col, ovsm_col_index.row) - - mat = mat.get(self.queue) - resample_mat = resampler.full_resample_matrix(self.queue).get(self.queue) + waa = _get_weights_and_area_elements(queue, source, source_where) + mat *= waa[index_set.linear_col_indices] - rmat = np.empty(ovsm_index_set.nblocks, dtype=np.object) - for i in range(ovsm_index_set.nblocks): - mat_blk = ovsm_col_index.block_take(mat, i) - resample_blk = ovsm_row_index.take(resample_mat, i) - - rmat[i] = mat_blk.dot(resample_blk).reshape(-1) - mat = np.hstack(rmat) - - # TODO:: multiply with rec_density + # TODO: multiply with rec_density return mat @@ -548,10 +532,19 @@ class FarFieldBlockBuilder(MatrixBlockBuilderBase): return np.equal(tgtindices, srcindices).astype(np.float64) def map_int_g(self, expr): - source = self.places[expr.source] - target_discr = self.places[expr.target] + where_source = expr.source + if where_source is DEFAULT_SOURCE: + where_source = sym.QBXSourceQuadStage2(where_source) + + where_target = expr.target + if where_target is DEFAULT_TARGET + where_target = sym.QBXSourceQuadStage2(expr.target) + + from pytential.symbolic.execution import _get_discretization + source, source_discr = _get_discretization(self.places, where_source) + _, target_discr = _get_discretization(self.places, where_target) - if source.density_discr is not target_discr: + if source_discr is not target_discr: raise NotImplementedError() rec_density = self.dummy.rec(expr.density) @@ -562,11 +555,8 @@ class FarFieldBlockBuilder(MatrixBlockBuilderBase): if len(rec_density.shape) != 2: raise NotImplementedError("layer potentials on non-variables") - try: - kernel = expr.kernel.inner_kernel - except AttributeError: - kernel = expr.kernel - kernel_args = _get_kernel_args(self, kernel, expr, source) + kernel_args = _get_kernel_args(self, + kernel.get_base_kernel(), expr, source) if self.exclude_self: kernel_args["target_to_source"] = \ cl.array.arange(self.queue, 0, target_discr.nnodes, dtype=np.int) @@ -577,7 +567,7 @@ class FarFieldBlockBuilder(MatrixBlockBuilderBase): _, (mat,) = mat_gen(self.queue, targets=target_discr.nodes(), - sources=source.density_discr.nodes(), + sources=source_discr.nodes(), index_set=self.index_set, **kernel_args) mat = mat.get() diff --git a/test/test_matrix.py b/test/test_matrix.py index b266c2e8..f36614fd 100644 --- a/test/test_matrix.py +++ b/test/test_matrix.py @@ -330,8 +330,7 @@ def test_qbx_block_builder(ctx_factory, ndim, lpot_id, visualize=False): @pytest.mark.parametrize('where', - [None, - (DEFAULT_SOURCE, DEFAULT_TARGET), + [(DEFAULT_SOURCE, DEFAULT_TARGET), (QBXSourceStage1(DEFAULT_SOURCE), QBXSourceStage1(DEFAULT_TARGET)), (QBXSourceQuadStage2(DEFAULT_SOURCE), @@ -391,22 +390,12 @@ def test_build_matrix_where(ctx_factory, where, visualize=False): from pytential.symbolic.execution import build_matrix mat = build_matrix(queue, qbx, op, u_sym, auto_where=where) - if where is None: - source_where, target_where = DEFAULT_SOURCE, DEFAULT_TARGET - else: - source_where, target_where = where - - if isinstance(source_where, QBXSourceQuadStage2): - n = qbx.quad_stage2_density_discr.nnodes - else: - n = qbx.density_discr.nnodes - - if isinstance(target_where, QBXSourceQuadStage2): - m = qbx.quad_stage2_density_discr.nnodes - else: - m = qbx.density_discr.nnodes + from pytential.symbolic.execution import _get_discretization + source_where, target_where = where + _, source_discr = _get_discretization(qbx, source_where) + _, target_discr = _get_discretization(qbx, target_where) - assert mat.shape == (m, n) + assert mat.shape == (target_discr.nnodes, source_discr.nnodes) if __name__ == "__main__": -- GitLab From c7dceec90ce02416cef542cdd3f22a1fc3f29e62 Mon Sep 17 00:00:00 2001 From: Alex Fikl Date: Tue, 31 Jul 2018 16:17:27 -0500 Subject: [PATCH 13/34] matrix: fix some tests and bugs --- pytential/symbolic/matrix.py | 54 ++++++++++++++++++++---------------- test/test_matrix.py | 13 +++++---- 2 files changed, 38 insertions(+), 29 deletions(-) diff --git a/pytential/symbolic/matrix.py b/pytential/symbolic/matrix.py index 4e6ded89..80f23296 100644 --- a/pytential/symbolic/matrix.py +++ b/pytential/symbolic/matrix.py @@ -84,7 +84,7 @@ def _get_kernel_args(mapper, kernel, expr, source): return kernel_args -def _weights_and_area_elements(queue, source, where): +def _get_weights_and_area_elements(queue, source, where): if isinstance(where, sym.QBXSourceQuadStage2): return source.weights_and_area_elements().with_queue(queue) @@ -105,12 +105,17 @@ def _get_centers_and_expansion_radii(queue, source, target_discr, qbx_forced_lim radii = source._expansion_radii('nsources') from pytential.qbx.target_assoc import associate_targets_to_qbx_centers + if source.density_discr is target_discr: + target_association_tolerance = 1.0e-10 + else: + target_association_tolerance = 1.0e-1 + code_container = source.target_association_code_container assoc = associate_targets_to_qbx_centers( source, code_container.get_wrangler(queue), [(target_discr, qbx_forced_limit)], - target_association_tolerance=1.0e-1) + target_association_tolerance=target_association_tolerance) centers = [cl.array.take(c, assoc.target_to_center, queue=queue) for c in centers] @@ -226,11 +231,11 @@ class MatrixBuilder(EvaluationMapperBase): def map_int_g(self, expr): where_source = expr.source - if where_source is DEFAULT_SOURCE: + if where_source is sym.DEFAULT_SOURCE: where_source = sym.QBXSourceQuadStage2(where_source) where_target = expr.target - if where_target is DEFAULT_TARGET + if where_target is sym.DEFAULT_TARGET: where_target = sym.QBXSourceStage1(expr.target) from pytential.symbolic.execution import _get_discretization @@ -268,7 +273,7 @@ class MatrixBuilder(EvaluationMapperBase): **kernel_args) mat = mat.get() - waa = _weights_and_area_elements(self.queue, source, where_source) + waa = _get_weights_and_area_elements(self.queue, source, where_source) mat[:, :] *= waa.get(self.queue) if target_discr.nnodes != source_discr.nnodes: @@ -332,11 +337,11 @@ class P2PMatrixBuilder(MatrixBuilder): def map_int_g(self, expr): where_source = expr.source - if where_source is DEFAULT_SOURCE: - where_source = sym.QBXSourceQuadStage2(where_source) + if where_source is sym.DEFAULT_SOURCE: + where_source = sym.QBXSourceStage1(where_source) where_target = expr.target - if where_target is DEFAULT_TARGET + if where_target is sym.DEFAULT_TARGET: where_target = sym.QBXSourceStage1(expr.target) from pytential.symbolic.execution import _get_discretization @@ -352,8 +357,8 @@ class P2PMatrixBuilder(MatrixBuilder): if len(rec_density.shape) != 2: raise NotImplementedError("layer potentials on non-variables") - kernel_args = _get_kernel_args(self, - kernel.get_base_kernel(), expr, source) + kernel = expr.kernel.get_base_kernel() + kernel_args = _get_kernel_args(self, kernel, expr, source) if self.exclude_self: kernel_args["target_to_source"] = \ cl.array.arange(self.queue, 0, target_discr.nnodes, dtype=np.int) @@ -462,12 +467,12 @@ class NearFieldBlockBuilder(MatrixBlockBuilderBase): def map_int_g(self, expr): where_source = expr.source - if where_source is DEFAULT_SOURCE: - where_source = sym.QBXSourceQuadStage2(where_source) + if where_source is sym.DEFAULT_SOURCE: + where_source = sym.QBXSourceStage1(where_source) where_target = expr.target - if where_target is DEFAULT_TARGET - where_target = sym.QBXSourceQuadStage2(expr.target) + if where_target is sym.DEFAULT_TARGET: + where_target = sym.QBXSourceStage1(expr.target) from pytential.symbolic.execution import _get_discretization source, source_discr = _get_discretization(self.places, where_source) @@ -496,7 +501,7 @@ class NearFieldBlockBuilder(MatrixBlockBuilderBase): from pytential.qbx.utils import get_centers_on_side assert abs(expr.qbx_forced_limit) > 0 - centers, radii = _get_centers_and_expansion_radii(queue, + centers, radii = _get_centers_and_expansion_radii(self.queue, source, target_discr, expr.qbx_forced_limit) _, (mat,) = mat_gen(self.queue, @@ -504,11 +509,12 @@ class NearFieldBlockBuilder(MatrixBlockBuilderBase): sources=source_discr.nodes(), centers=centers, expansion_radii=radii, - index_set=index_set, + index_set=self.index_set, **kernel_args) - waa = _get_weights_and_area_elements(queue, source, source_where) - mat *= waa[index_set.linear_col_indices] + waa = _get_weights_and_area_elements(self.queue, source, where_source) + mat *= waa[self.index_set.linear_col_indices] + mat = mat.get(self.queue) # TODO: multiply with rec_density @@ -533,12 +539,12 @@ class FarFieldBlockBuilder(MatrixBlockBuilderBase): def map_int_g(self, expr): where_source = expr.source - if where_source is DEFAULT_SOURCE: - where_source = sym.QBXSourceQuadStage2(where_source) + if where_source is sym.DEFAULT_SOURCE: + where_source = sym.QBXSourceStage1(where_source) where_target = expr.target - if where_target is DEFAULT_TARGET - where_target = sym.QBXSourceQuadStage2(expr.target) + if where_target is sym.DEFAULT_TARGET: + where_target = sym.QBXSourceStage1(expr.target) from pytential.symbolic.execution import _get_discretization source, source_discr = _get_discretization(self.places, where_source) @@ -555,8 +561,8 @@ class FarFieldBlockBuilder(MatrixBlockBuilderBase): if len(rec_density.shape) != 2: raise NotImplementedError("layer potentials on non-variables") - kernel_args = _get_kernel_args(self, - kernel.get_base_kernel(), expr, source) + kernel = expr.kernel.get_base_kernel() + kernel_args = _get_kernel_args(self, kernel, expr, source) if self.exclude_self: kernel_args["target_to_source"] = \ cl.array.arange(self.queue, 0, target_discr.nnodes, dtype=np.int) diff --git a/test/test_matrix.py b/test/test_matrix.py index f36614fd..1629ab5d 100644 --- a/test/test_matrix.py +++ b/test/test_matrix.py @@ -273,8 +273,10 @@ def test_qbx_block_builder(ctx_factory, ndim, lpot_id, visualize=False): srcindices = _build_block_index(qbx.density_discr) from pytential.symbolic.execution import prepare_places, prepare_expr - places = prepare_places(qbx) - expr = prepare_expr(places, op) + where = (QBXSourceStage1(DEFAULT_SOURCE), + QBXSourceStage1(DEFAULT_TARGET)) + places = prepare_places(qbx, auto_where=where) + expr = prepare_expr(places, op, auto_where=where) from sumpy.tools import MatrixBlockIndexRanges index_set = MatrixBlockIndexRanges(ctx, tgtindices, srcindices) @@ -283,18 +285,19 @@ def test_qbx_block_builder(ctx_factory, ndim, lpot_id, visualize=False): mbuilder = NearFieldBlockBuilder(queue, dep_expr=u_sym, other_dep_exprs=[], - dep_source=places[DEFAULT_SOURCE], + dep_source=places[where[0]], places=places, context={}, index_set=index_set) blk = mbuilder(expr) + from pytential.symbolic.execution import _get_discretization from pytential.symbolic.matrix import MatrixBuilder mbuilder = MatrixBuilder(queue, dep_expr=u_sym, other_dep_exprs=[], - dep_source=places[DEFAULT_SOURCE], - dep_discr=places[DEFAULT_SOURCE].density_discr, + dep_source=places[where[0]], + dep_discr=_get_discretization(places, where[0])[1], places=places, context={}) mat = mbuilder(expr) -- GitLab From 5c541f642e9fcccc186eecee924473679342ae5d Mon Sep 17 00:00:00 2001 From: Alex Fikl Date: Tue, 31 Jul 2018 18:34:17 -0500 Subject: [PATCH 14/34] matrix: flake8 fixes --- pytential/symbolic/matrix.py | 3 --- test/test_matrix.py | 4 ++-- 2 files changed, 2 insertions(+), 5 deletions(-) diff --git a/pytential/symbolic/matrix.py b/pytential/symbolic/matrix.py index 80f23296..200e1e99 100644 --- a/pytential/symbolic/matrix.py +++ b/pytential/symbolic/matrix.py @@ -29,8 +29,6 @@ import pyopencl.array # noqa import six from six.moves import intern -from pytools import memoize_method - from pytential.symbolic.mappers import EvaluationMapperBase import pytential.symbolic.primitives as sym from pytential.symbolic.execution import bind @@ -499,7 +497,6 @@ class NearFieldBlockBuilder(MatrixBlockBuilderBase): mat_gen = LayerPotentialMatrixBlockGenerator( self.queue.context, (local_expn,)) - from pytential.qbx.utils import get_centers_on_side assert abs(expr.qbx_forced_limit) > 0 centers, radii = _get_centers_and_expansion_radii(self.queue, source, target_discr, expr.qbx_forced_limit) diff --git a/test/test_matrix.py b/test/test_matrix.py index 1629ab5d..842029b1 100644 --- a/test/test_matrix.py +++ b/test/test_matrix.py @@ -379,8 +379,8 @@ def test_build_matrix_where(ctx_factory, where, visualize=False): continue pt.figure(figsize=(10, 8), dpi=300) - pt.plot(sources[0], sources[1], 'k'); - pt.plot(targets[0], targets[1], 'ko'); + pt.plot(sources[0], sources[1], 'k') + pt.plot(targets[0], targets[1], 'ko') line = pt.plot(targets[0, itgt], targets[1, itgt], 'o') c = pt.Circle([centers[0][i], centers[1][i]], radii[i // 2], -- GitLab From a39fda90acd53cdab377baf8620ddf149de93f68 Mon Sep 17 00:00:00 2001 From: Alex Fikl Date: Tue, 31 Jul 2018 19:52:28 -0500 Subject: [PATCH 15/34] execution: add some docs --- pytential/symbolic/execution.py | 77 ++++++++++++++++++++++++++++++--- 1 file changed, 72 insertions(+), 5 deletions(-) diff --git a/pytential/symbolic/execution.py b/pytential/symbolic/execution.py index 2ec5b05f..17a09fd7 100644 --- a/pytential/symbolic/execution.py +++ b/pytential/symbolic/execution.py @@ -284,22 +284,42 @@ class MatVecOp: # {{{ expression prep def _domains_default(nresults, places, domains, default_val): + """ + :arg nresults: number of results. + :arg places: result of :func:`prepare_places`. + :arg domains: recommended domains. + :arg default_val: default value for domains which are not provided. + + :return: a list of domains for each result. If domains is `None`, each + element in the list is *default_val*. If *domains* is not a list + of domains, each element in the resulting list is *domains*. Otherwise, + *domains* is returned as is. + """ if domains is None: if default_val not in places: raise RuntimeError("'domains is None' requires " "default domain to be defined") dom_name = default_val return nresults * [dom_name] - elif not isinstance(domains, (list, tuple)): dom_name = domains return nresults * [dom_name] - else: + assert len(domains) == nresults return domains def _where_default(places, auto_where=None): + """ + :arg places: result of :func:`prepare_places`. + :arg auto_where: identifiers for source and/or target locations. If `None`, + `where` attributes are automatically found. + + :return: if *auto_where* is provided, it is returned as is. If *places* + was obtained from :func:`prepare_places`, the default is given by its + keys. Otherwise, a tuple of `(DEFAULT_SOURCE, DEFAULT_TARGET)` is + returned. + """ if auto_where is None: if not isinstance(places, dict): return DEFAULT_SOURCE, DEFAULT_TARGET @@ -311,6 +331,18 @@ def _where_default(places, auto_where=None): def prepare_places(places, auto_where=None): + """ + :arg places: a mapping of symbolic names to + :class:`~meshmode.discretization.Discretization` objects or a subclass + of :class:`~pytential.target.TargetBase`. + :arg auto_where: identifiers for source and/or target locations. If `None`, + `where` attributes are automatically found. + + :return: a mapping of symbolic names, same as the input if it was already + such a mapping. If not, a mapping is constructed using the values + of *auto_where*, if provided, as keys and appropriate discretization + objects as values. + """ from meshmode.discretization import Discretization from pytential.source import LayerPotentialSourceBase from pytential.target import TargetBase @@ -351,9 +383,13 @@ def prepare_places(places, auto_where=None): def prepare_expr(places, expr, auto_where=None): """ - :arg places: result of :func:`prepare_places` - :arg auto_where: For simple source-to-self or source-to-target - evaluations, find 'where' attributes automatically. + :arg places: result of :func:`prepare_places`. + :arg expr: an array of symbolic expressions. + :arg auto_where: identifiers for source and/or target locations. If `None`, + `where` attributes are automatically found. + + :return: processed symbolic expressions, tagger with the given `where` + identifiers. """ from pytential.source import LayerPotentialSourceBase @@ -375,6 +411,25 @@ def prepare_expr(places, expr, auto_where=None): def prepare_expression(places, exprs, input_exprs, domains=None, auto_where=None): + """ + :arg places: a mapping of symbolic names to + :class:`~meshmode.discretization.Discretization` objects or a subclass + of :class:`~pytential.target.TargetBase`. + :arg exprs: an array or a single symbolic expression. + :arg input_exprs: an array or a single symbolic expression that is taken + as input by *exprs*. + :arg domains: a list of discretization identifiers, indicating the domains + on which the inputs live. If given, each element of the list must be + a key in mapping *places* and correspond to an *auto_where* + identifier. + :arg auto_where: identifiers for source and/or target locations. If `None`, + `where` attributes are automatically found. + + :return: a tuple of `(places, exprs, input_exprs, domains)`, where each + element was appropriately processed so that it can be used by + :class:`BoundExpression`, :func:`build_matrix`, etc. + """ + auto_where = _where_default(places, auto_where) places = prepare_places(places, auto_where=auto_where) exprs = prepare_expr(places, exprs, auto_where=auto_where) @@ -398,6 +453,18 @@ def prepare_expression(places, exprs, input_exprs, # {{{ bound expression def _get_discretization(places, where, default_source=QBXSourceStage1): + """ + :arg places: a mapping of symbolic names to + :class:`~meshmode.discretization.Discretization` objects or a subclass + of :class:`~pytential.target.TargetBase`. + :arg where: identifier for source or target locations. + :arg default_source: specific source location in case `where` is + :class:`pytential.symbolic.primitives.DEFAULT_SOURCE` or + :class:`pytential.symbolic.primitives.DEFAULT_TARGET`. + + :return: a :class:`~meshmode.discretization.Discretization`, from + *places* corresponding to *where*. + """ from pytential.source import LayerPotentialSourceBase if where is DEFAULT_SOURCE or where is DEFAULT_TARGET: -- GitLab From a2d1c9773a719d8dea80cea24fdc0bb98e1d087e Mon Sep 17 00:00:00 2001 From: Alex Fikl Date: Fri, 3 Aug 2018 19:44:14 -0500 Subject: [PATCH 16/34] execution: first try at wrapping places handling in a class --- pytential/symbolic/execution.py | 64 ++++++++++++++++++++++++++++++--- 1 file changed, 60 insertions(+), 4 deletions(-) diff --git a/pytential/symbolic/execution.py b/pytential/symbolic/execution.py index 17a09fd7..d9f64ee8 100644 --- a/pytential/symbolic/execution.py +++ b/pytential/symbolic/execution.py @@ -452,6 +452,64 @@ def prepare_expression(places, exprs, input_exprs, # {{{ bound expression +class BindingLocation(object): + def __init__(self, places, auto_where=None): + from meshmode.discretization import Discretization + from pytential.source import LayerPotentialSourceBase, PotentialSource + from pytential.target import TargetBase + + if auto_where is None: + auto_where = DEFAULT_SOURCE, DEFAULT_TARGET + where_source, where_target = auto_where + + def get_lpot_discr(where): + if where is DEFAULT_SOURCE or where is DEFAULT_TARGET: + where = QBXSourceStage1(where) + + if isinstance(where, QBXSourceStage1): + return places.density_discr + if isinstance(where, QBXSourceStage2): + return places.stage2_density_discr + if isinstance(where, QBXSourceQuadStage2): + return places.quad_stage2_density_discr + + return None + + self.places = {} + if isinstance(places, LayerPotentialSourceBase): + self[where_source] = get_lpot_discr(places, where_source) + self[where_target] = get_lpot_discr(places, where_target) + elif isinstance(places, (Discretization, TargetBase)): + self[where_target] = places + elif isinstance(places, tuple): + source_discr, target_discr = places + self[where_source] = source_discr + self[where_target] = target_discr + else: + for key, value in six.iteritems(places): + self[key] = value + + for p in six.itervalues(self.places): + if not isinstance(p, (PotentialSource, TargetBase, Discretization)): + raise TypeError("Must pass discretization, targets or " + "layer potential sources as 'places'") + + self.source = None + if isinstance(places, LayerPotentialSourceBase): + self.source = places + + self.caches = {} + + def __getitem__(self, where): + return self.places[where] + + def __setitem__(self, where, discr): + self.places[where] = discr + + def get_cache(self, name): + return self.caches.setdefault(name, {}) + + def _get_discretization(places, where, default_source=QBXSourceStage1): """ :arg places: a mapping of symbolic names to @@ -491,12 +549,10 @@ def _get_discretization(places, where, default_source=QBXSourceStage1): return lpot, discr -class BoundExpression: - def __init__(self, places, sym_op_expr, sym_op_args=None): +class BoundExpression(object): + def __init__(self, places, sym_op_expr): self.places = places self.sym_op_expr = sym_op_expr - self.sym_op_args = sym_op_args - self.caches = {} from pytential.symbolic.compiler import OperatorCompiler -- GitLab From 04092e6904d1ee6a882db52eb871db4ae856466e Mon Sep 17 00:00:00 2001 From: Alex Fikl Date: Fri, 3 Aug 2018 19:50:35 -0500 Subject: [PATCH 17/34] matrix: add a note about target association tolerance --- pytential/symbolic/matrix.py | 62 ++++++++++++++++++++---------------- 1 file changed, 35 insertions(+), 27 deletions(-) diff --git a/pytential/symbolic/matrix.py b/pytential/symbolic/matrix.py index 200e1e99..63691696 100644 --- a/pytential/symbolic/matrix.py +++ b/pytential/symbolic/matrix.py @@ -84,42 +84,50 @@ def _get_kernel_args(mapper, kernel, expr, source): def _get_weights_and_area_elements(queue, source, where): if isinstance(where, sym.QBXSourceQuadStage2): - return source.weights_and_area_elements().with_queue(queue) - - # NOTE: copied from `weights_and_area_elements`, but using the - # discretization given by `where` and no interpolation - from pytential.symbolic.execution import _get_discretization - discr = _get_discretization(source, where) + waa = source.weights_and_area_elements().with_queue(queue) + else: + # NOTE: copied from `weights_and_area_elements`, but using the + # discretization given by `where` and no interpolation + from pytential.symbolic.execution import _get_discretization + discr = _get_discretization(source, where) - area = bind(discr, sym.area_element(source.ambient_dim, source.dim))(queue) - qweight = bind(discr, sym.QWeight())(queue) + area = bind(discr, sym.area_element(source.ambient_dim, source.dim))(queue) + qweight = bind(discr, sym.QWeight())(queue) + waa = area * qweight - return area * qweight + return waa def _get_centers_and_expansion_radii(queue, source, target_discr, qbx_forced_limit): - from pytential.qbx.utils import get_interleaved_centers - centers = get_interleaved_centers(queue, source) - radii = source._expansion_radii('nsources') - - from pytential.qbx.target_assoc import associate_targets_to_qbx_centers + # NOTE: skip expensive target association if source.density_discr is target_discr: - target_association_tolerance = 1.0e-10 + from pytential.qbx.utils import get_centers_on_side + centers = get_centers_on_side(source, qbx_forced_limit) + radii = source._expansion_radii('nsources') else: + from pytential.qbx.utils import get_interleaved_centers + centers = get_interleaved_centers(queue, source) + radii = source._expansion_radii('nsources') + + # NOTE: using a very small tolerance to make sure all the stage2 + # targets are associated to a center. We can't use the user provided + # source.target_association_tolerance here because it will likely be + # way too small. target_association_tolerance = 1.0e-1 - code_container = source.target_association_code_container - assoc = associate_targets_to_qbx_centers( - source, - code_container.get_wrangler(queue), - [(target_discr, qbx_forced_limit)], - target_association_tolerance=target_association_tolerance) - - centers = [cl.array.take(c, assoc.target_to_center, queue=queue) - for c in centers] - radii = cl.array.take(radii, - (assoc.target_to_center.with_queue(queue) / 2.0).astype(np.int), - queue=queue) + from pytential.qbx.target_assoc import associate_targets_to_qbx_centers + code_container = source.target_association_code_container + assoc = associate_targets_to_qbx_centers( + source, + code_container.get_wrangler(queue), + [(target_discr, qbx_forced_limit)], + target_association_tolerance=target_association_tolerance) + + centers = [cl.array.take(c, assoc.target_to_center, queue=queue) + for c in centers] + radii = cl.array.take(radii, + (assoc.target_to_center.with_queue(queue) / 2.0).astype(np.int), + queue=queue) return centers, radii -- GitLab From 05149897384bb16df446ac91f5781037556b5281 Mon Sep 17 00:00:00 2001 From: Alex Fikl Date: Sun, 5 Aug 2018 20:02:01 -0500 Subject: [PATCH 18/34] execution: use the new places class everywhere --- pytential/symbolic/compiler.py | 4 +- pytential/symbolic/execution.py | 291 ++++++++++---------------------- pytential/symbolic/mappers.py | 5 +- pytential/symbolic/matrix.py | 124 +++++--------- test/test_matrix.py | 63 ++++--- 5 files changed, 175 insertions(+), 312 deletions(-) diff --git a/pytential/symbolic/compiler.py b/pytential/symbolic/compiler.py index 17d23ebf..ecb1fdca 100644 --- a/pytential/symbolic/compiler.py +++ b/pytential/symbolic/compiler.py @@ -224,7 +224,7 @@ class ComputePotentialInstruction(Instruction): ", ".join(args), "\n ".join(lines)) def get_exec_function(self, exec_mapper): - source = exec_mapper.bound_expr.places[self.source] + source = exec_mapper.bound_expr.places.lpot return source.exec_compute_potential_insn def __hash__(self): @@ -444,7 +444,7 @@ class OperatorCompiler(IdentityMapper): def op_group_features(self, expr): from pytential.symbolic.primitives import hashable_kernel_args return ( - self.places[expr.source].op_group_features(expr) + self.places.lpot.op_group_features(expr) + hashable_kernel_args(expr.kernel_arguments)) @memoize_method diff --git a/pytential/symbolic/execution.py b/pytential/symbolic/execution.py index d9f64ee8..0217be3d 100644 --- a/pytential/symbolic/execution.py +++ b/pytential/symbolic/execution.py @@ -283,7 +283,7 @@ class MatVecOp: # {{{ expression prep -def _domains_default(nresults, places, domains, default_val): +def _prepare_domains(nresults, places, domains, default_val): """ :arg nresults: number of results. :arg places: result of :func:`prepare_places`. @@ -309,79 +309,7 @@ def _domains_default(nresults, places, domains, default_val): return domains -def _where_default(places, auto_where=None): - """ - :arg places: result of :func:`prepare_places`. - :arg auto_where: identifiers for source and/or target locations. If `None`, - `where` attributes are automatically found. - - :return: if *auto_where* is provided, it is returned as is. If *places* - was obtained from :func:`prepare_places`, the default is given by its - keys. Otherwise, a tuple of `(DEFAULT_SOURCE, DEFAULT_TARGET)` is - returned. - """ - if auto_where is None: - if not isinstance(places, dict): - return DEFAULT_SOURCE, DEFAULT_TARGET - if DEFAULT_TARGET in places: - return DEFAULT_SOURCE, DEFAULT_TARGET - return tuple(places.keys()) - - return auto_where - - -def prepare_places(places, auto_where=None): - """ - :arg places: a mapping of symbolic names to - :class:`~meshmode.discretization.Discretization` objects or a subclass - of :class:`~pytential.target.TargetBase`. - :arg auto_where: identifiers for source and/or target locations. If `None`, - `where` attributes are automatically found. - - :return: a mapping of symbolic names, same as the input if it was already - such a mapping. If not, a mapping is constructed using the values - of *auto_where*, if provided, as keys and appropriate discretization - objects as values. - """ - from meshmode.discretization import Discretization - from pytential.source import LayerPotentialSourceBase - from pytential.target import TargetBase - - where_source, where_target = _where_default(places, auto_where=auto_where) - if isinstance(places, LayerPotentialSourceBase): - _, target_discr = _get_discretization(places, where_target) - places = { - where_source: places, - where_target: target_discr, - } - elif isinstance(places, (Discretization, TargetBase)): - places = { - where_target: places, - } - - elif isinstance(places, tuple): - source_discr, target_discr = places - places = { - where_source: source_discr, - where_target: target_discr, - } - del source_discr - del target_discr - - def cast_to_place(discr): - from pytential.target import TargetBase - from pytential.source import PotentialSource - if not isinstance(discr, (Discretization, TargetBase, PotentialSource)): - raise TypeError("must pass discretizations, " - "layer potential sources or targets as 'places'") - return discr - - return dict( - (key, cast_to_place(value)) - for key, value in six.iteritems(places)) - - -def prepare_expr(places, expr, auto_where=None): +def _prepare_expr(places, expr): """ :arg places: result of :func:`prepare_places`. :arg expr: an array of symbolic expressions. @@ -396,57 +324,13 @@ def prepare_expr(places, expr, auto_where=None): from pytential.symbolic.mappers import ( ToTargetTagger, DerivativeBinder) - auto_where = _where_default(places, auto_where=auto_where) - if auto_where: - expr = ToTargetTagger(*auto_where)(expr) - + expr = ToTargetTagger(*places.where)(expr) expr = DerivativeBinder()(expr) - - for name, place in six.iteritems(places): - if isinstance(place, LayerPotentialSourceBase): - expr = place.preprocess_optemplate(name, places, expr) + if isinstance(places.lpot, LayerPotentialSourceBase): + expr = places.lpot.preprocess_optemplate(places.source, places, expr) return expr - -def prepare_expression(places, exprs, input_exprs, - domains=None, auto_where=None): - """ - :arg places: a mapping of symbolic names to - :class:`~meshmode.discretization.Discretization` objects or a subclass - of :class:`~pytential.target.TargetBase`. - :arg exprs: an array or a single symbolic expression. - :arg input_exprs: an array or a single symbolic expression that is taken - as input by *exprs*. - :arg domains: a list of discretization identifiers, indicating the domains - on which the inputs live. If given, each element of the list must be - a key in mapping *places* and correspond to an *auto_where* - identifier. - :arg auto_where: identifiers for source and/or target locations. If `None`, - `where` attributes are automatically found. - - :return: a tuple of `(places, exprs, input_exprs, domains)`, where each - element was appropriately processed so that it can be used by - :class:`BoundExpression`, :func:`build_matrix`, etc. - """ - - auto_where = _where_default(places, auto_where) - places = prepare_places(places, auto_where=auto_where) - exprs = prepare_expr(places, exprs, auto_where=auto_where) - - from pytools.obj_array import is_obj_array, make_obj_array - if not is_obj_array(exprs): - exprs = make_obj_array([exprs]) - try: - input_exprs = list(input_exprs) - except TypeError: - # not iterable, wrap in a list - input_exprs = [input_exprs] - - domains = _domains_default(len(input_exprs), places, domains, auto_where[0]) - - return places, exprs, input_exprs, domains - # }}} @@ -460,93 +344,91 @@ class BindingLocation(object): if auto_where is None: auto_where = DEFAULT_SOURCE, DEFAULT_TARGET - where_source, where_target = auto_where - - def get_lpot_discr(where): - if where is DEFAULT_SOURCE or where is DEFAULT_TARGET: - where = QBXSourceStage1(where) - - if isinstance(where, QBXSourceStage1): - return places.density_discr - if isinstance(where, QBXSourceStage2): - return places.stage2_density_discr - if isinstance(where, QBXSourceQuadStage2): - return places.quad_stage2_density_discr - - return None + self.where = auto_where + self.lpot = None self.places = {} if isinstance(places, LayerPotentialSourceBase): - self[where_source] = get_lpot_discr(places, where_source) - self[where_target] = get_lpot_discr(places, where_target) + self.lpot = places + self.places[self.source] = self[self.source] + self.places[self.target] = self[self.target] elif isinstance(places, (Discretization, TargetBase)): - self[where_target] = places + self.places[self.target] = places elif isinstance(places, tuple): source_discr, target_discr = places - self[where_source] = source_discr - self[where_target] = target_discr + if isinstance(source_discr, PotentialSource): + self.lpot = source_discr + if isinstance(source_discr, LayerPotentialSourceBase): + source_discr = self[self.source] + + self.places[self.source] = source_discr + self.places[self.target] = target_discr else: for key, value in six.iteritems(places): - self[key] = value + self.places[key] = value for p in six.itervalues(self.places): if not isinstance(p, (PotentialSource, TargetBase, Discretization)): raise TypeError("Must pass discretization, targets or " "layer potential sources as 'places'") - self.source = None - if isinstance(places, LayerPotentialSourceBase): - self.source = places - self.caches = {} + @property + def source(self): + return self.where[0] + + @property + def target(self): + return self.where[1] + + def _get_lpot_discretization(self, where): + from pytential.source import LayerPotentialSourceBase + if not isinstance(self.lpot, LayerPotentialSourceBase): + return None + + if where is DEFAULT_SOURCE or where is DEFAULT_TARGET: + where = QBXSourceStage1(where) + + if isinstance(where, QBXSourceStage1): + return self.lpot.density_discr + if isinstance(where, QBXSourceStage2): + return self.lpot.stage2_density_discr + if isinstance(where, QBXSourceQuadStage2): + return self.lpot.quad_stage2_density_discr + + return None + def __getitem__(self, where): - return self.places[where] + if where in self.places: + discr = self.places[where] + else: + discr = self._get_lpot_discretization(where) - def __setitem__(self, where, discr): - self.places[where] = discr + if discr is None: + from pytential.symbolic.mappers import stringify_where + raise KeyError("Identifier '{}' not in places".format( + stringify_where(where))) - def get_cache(self, name): - return self.caches.setdefault(name, {}) + return discr + def __iter__(self): + return iter(self.places.keys()) -def _get_discretization(places, where, default_source=QBXSourceStage1): - """ - :arg places: a mapping of symbolic names to - :class:`~meshmode.discretization.Discretization` objects or a subclass - of :class:`~pytential.target.TargetBase`. - :arg where: identifier for source or target locations. - :arg default_source: specific source location in case `where` is - :class:`pytential.symbolic.primitives.DEFAULT_SOURCE` or - :class:`pytential.symbolic.primitives.DEFAULT_TARGET`. - - :return: a :class:`~meshmode.discretization.Discretization`, from - *places* corresponding to *where*. - """ - from pytential.source import LayerPotentialSourceBase + def keys(self): + return self.places.keys() - if where is DEFAULT_SOURCE or where is DEFAULT_TARGET: - where = default_source(where) + def values(self): + return self.places.values() - if isinstance(places, LayerPotentialSourceBase): - lpot = places - else: - try: - lpot = places[where.where] - except KeyError: - lpot = places[where] - is_lpot = isinstance(lpot, LayerPotentialSourceBase) - - if isinstance(where, QBXSourceStage1): - discr = lpot.density_discr if is_lpot else lpot - elif isinstance(where, QBXSourceStage2): - discr = lpot.stage2_density_discr if is_lpot else lpot - elif isinstance(where, QBXSourceQuadStage2): - discr = lpot.quad_stage2_density_discr if is_lpot else lpot - else: - raise ValueError("Unknown 'where': {}".format(type(where))) + def items(self): + return self.places.items() - return lpot, discr + def __hash__(self): + return hash(tuple(self.values())) + + def get_cache(self, name): + return self.caches.setdefault(name, {}) class BoundExpression(object): @@ -559,11 +441,10 @@ class BoundExpression(object): self.code = OperatorCompiler(self.places)(sym_op_expr) def get_cache(self, name): - return self.caches.setdefault(name, {}) + return self.places.get_cache(name) def get_discretization(self, where): - _, discr = _get_discretization(self.places, where) - return discr + return self.places[where] def scipy_op(self, queue, arg_name, dtype, domains=None, **extra_args): """ @@ -581,7 +462,7 @@ class BoundExpression(object): else: nresults = 1 - domains = _domains_default(nresults, self.places, domains, + domains = _prepare_domains(nresults, self.places, domains, DEFAULT_TARGET) total_dofs = 0 @@ -617,8 +498,9 @@ def bind(places, expr, auto_where=None): evaluations, find 'where' attributes automatically. """ - places = prepare_places(places) - expr = prepare_expr(places, expr, auto_where=auto_where) + places = BindingLocation(places, auto_where=auto_where) + expr = _prepare_expr(places, expr) + return BoundExpression(places, expr) # }}} @@ -648,34 +530,36 @@ def build_matrix(queue, places, exprs, input_exprs, domains=None, evaluations, find 'where' attributes automatically. """ - from pytools import single_valued - from pytential.symbolic.matrix import MatrixBuilder, is_zero - if context is None: context = {} - auto_where = _where_default(places, auto_where=auto_where) - places, exprs, input_exprs, domains = \ - prepare_expression(places, exprs, input_exprs, - domains=domains, - auto_where=auto_where) + from pytools.obj_array import is_obj_array, make_obj_array + places = BindingLocation(places, auto_where=auto_where) + exprs = _prepare_expr(places, exprs) + if not is_obj_array(exprs): + exprs = make_obj_array([exprs]) + try: + input_exprs = list(input_exprs) + except TypeError: + # not iterable, wrap in a list + input_exprs = [input_exprs] + + domains = _prepare_domains(len(input_exprs), places, domains, places.source) + + from pytential.symbolic.matrix import MatrixBuilder, is_zero nblock_rows = len(exprs) nblock_columns = len(input_exprs) blocks = np.zeros((nblock_rows, nblock_columns), dtype=np.object) dtypes = [] for ibcol in range(nblock_columns): - dep_source, dep_discr = \ - _get_discretization(places, domains[ibcol]) - mbuilder = MatrixBuilder( queue, dep_expr=input_exprs[ibcol], other_dep_exprs=(input_exprs[:ibcol] + input_exprs[ibcol + 1:]), - dep_source=dep_source, - dep_discr=dep_discr, + dep_discr=places[domains[ibcol]], places=places, context=context) @@ -687,6 +571,7 @@ def build_matrix(queue, places, exprs, input_exprs, domains=None, if isinstance(block, np.ndarray): dtypes.append(block.dtype) + from pytools import single_valued block_row_counts = [ single_valued( blocks[ibrow, ibcol].shape[0] diff --git a/pytential/symbolic/mappers.py b/pytential/symbolic/mappers.py index 12b87894..a8696efe 100644 --- a/pytential/symbolic/mappers.py +++ b/pytential/symbolic/mappers.py @@ -399,9 +399,8 @@ class QBXPreprocessor(IdentityMapper): self.places = places def map_int_g(self, expr): - from pytential.symbolic.execution import _get_discretization - source, source_discr = _get_discretization(self.places, self.source_name) - _, target_discr = _get_discretization(self.places, expr.target) + source_discr = self.places[self.source_name] + target_discr = self.places[expr.target] if expr.qbx_forced_limit == 0: raise ValueError("qbx_forced_limit == 0 was a bad idea and " diff --git a/pytential/symbolic/matrix.py b/pytential/symbolic/matrix.py index 63691696..8c9a354e 100644 --- a/pytential/symbolic/matrix.py +++ b/pytential/symbolic/matrix.py @@ -82,25 +82,23 @@ def _get_kernel_args(mapper, kernel, expr, source): return kernel_args -def _get_weights_and_area_elements(queue, source, where): - if isinstance(where, sym.QBXSourceQuadStage2): +def _get_weights_and_area_elements(queue, source, source_discr): + if source.quad_stage2_density_discr is source_discr: waa = source.weights_and_area_elements().with_queue(queue) else: # NOTE: copied from `weights_and_area_elements`, but using the # discretization given by `where` and no interpolation - from pytential.symbolic.execution import _get_discretization - discr = _get_discretization(source, where) - - area = bind(discr, sym.area_element(source.ambient_dim, source.dim))(queue) - qweight = bind(discr, sym.QWeight())(queue) + area = bind(source_discr, + sym.area_element(source.ambient_dim, source.dim))(queue) + qweight = bind(source_discr, sym.QWeight())(queue) waa = area * qweight return waa def _get_centers_and_expansion_radii(queue, source, target_discr, qbx_forced_limit): - # NOTE: skip expensive target association if source.density_discr is target_discr: + # NOTE: skip expensive target association from pytential.qbx.utils import get_centers_on_side centers = get_centers_on_side(source, qbx_forced_limit) radii = source._expansion_radii('nsources') @@ -140,14 +138,13 @@ def _get_centers_and_expansion_radii(queue, source, target_discr, qbx_forced_lim # We'll cheat and build the matrix on the host. class MatrixBuilder(EvaluationMapperBase): - def __init__(self, queue, dep_expr, other_dep_exprs, dep_source, dep_discr, + def __init__(self, queue, dep_expr, other_dep_exprs, dep_discr, places, context): super(MatrixBuilder, self).__init__(context=context) self.queue = queue self.dep_expr = dep_expr self.other_dep_exprs = other_dep_exprs - self.dep_source = dep_source self.dep_discr = dep_discr self.places = places @@ -236,18 +233,13 @@ class MatrixBuilder(EvaluationMapperBase): return vecs_and_scalars def map_int_g(self, expr): + # TODO: should this go into QBXPreprocessor / LocationTagger where_source = expr.source if where_source is sym.DEFAULT_SOURCE: - where_source = sym.QBXSourceQuadStage2(where_source) - - where_target = expr.target - if where_target is sym.DEFAULT_TARGET: - where_target = sym.QBXSourceStage1(expr.target) + where_source = sym.QBXSourceQuadStage2(expr.source) - from pytential.symbolic.execution import _get_discretization - source, source_discr = _get_discretization(self.places, where_source) - _, target_discr = _get_discretization(self.places, where_target) - assert source_discr.nnodes >= target_discr.nnodes + source_discr = self.places[where_source] + target_discr = self.places[expr.target] rec_density = self.rec(expr.density) if is_zero(rec_density): @@ -258,10 +250,10 @@ class MatrixBuilder(EvaluationMapperBase): raise NotImplementedError("layer potentials on non-variables") kernel = expr.kernel - kernel_args = _get_layer_potential_args(self, expr, source) + kernel_args = _get_layer_potential_args(self, expr, self.places.lpot) from sumpy.expansion.local import LineTaylorLocalExpansion - local_expn = LineTaylorLocalExpansion(kernel, source.qbx_order) + local_expn = LineTaylorLocalExpansion(kernel, self.places.lpot.qbx_order) from sumpy.qbx import LayerPotentialMatrixGenerator mat_gen = LayerPotentialMatrixGenerator( @@ -269,7 +261,7 @@ class MatrixBuilder(EvaluationMapperBase): assert abs(expr.qbx_forced_limit) > 0 centers, radii = _get_centers_and_expansion_radii(self.queue, - source, target_discr, expr.qbx_forced_limit) + self.places.lpot, target_discr, expr.qbx_forced_limit) _, (mat,) = mat_gen(self.queue, targets=target_discr.nodes(), @@ -279,11 +271,13 @@ class MatrixBuilder(EvaluationMapperBase): **kernel_args) mat = mat.get() - waa = _get_weights_and_area_elements(self.queue, source, where_source) + waa = _get_weights_and_area_elements(self.queue, self.places.lpot, source_discr) mat[:, :] *= waa.get(self.queue) if target_discr.nnodes != source_discr.nnodes: - resampler = source.direct_resampler + assert target_discr.nnodes < source_discr.nnodes + + resampler = self.places.lpot.direct_resampler resample_mat = resampler.full_resample_matrix(self.queue).get(self.queue) mat = mat.dot(resample_mat) @@ -321,7 +315,7 @@ class MatrixBuilder(EvaluationMapperBase): rec_arg = cl.array.to_device(self.queue, rec_arg) op = expr.function(sym.var("u")) - result = bind(self.dep_source, op)(self.queue, u=rec_arg) + result = bind(self.places.lpot, op)(self.queue, u=rec_arg) if isinstance(result, cl.array.Array): result = result.get() @@ -334,26 +328,16 @@ class MatrixBuilder(EvaluationMapperBase): # {{{ p2p matrix builder class P2PMatrixBuilder(MatrixBuilder): - def __init__(self, queue, dep_expr, other_dep_exprs, dep_source, dep_discr, + def __init__(self, queue, dep_expr, other_dep_exprs, dep_discr, places, context, exclude_self=True): super(P2PMatrixBuilder, self).__init__(queue, - dep_expr, other_dep_exprs, dep_source, dep_discr, places, context) + dep_expr, other_dep_exprs, dep_discr, places, context) self.exclude_self = exclude_self def map_int_g(self, expr): - where_source = expr.source - if where_source is sym.DEFAULT_SOURCE: - where_source = sym.QBXSourceStage1(where_source) - - where_target = expr.target - if where_target is sym.DEFAULT_TARGET: - where_target = sym.QBXSourceStage1(expr.target) - - from pytential.symbolic.execution import _get_discretization - source, source_discr = _get_discretization(self.places, where_source) - _, target_discr = _get_discretization(self.places, where_target) - assert source_discr.nnodes >= target_discr.nnodes + source_discr = self.places[expr.source] + target_discr = self.places[expr.target] rec_density = self.rec(expr.density) if is_zero(rec_density): @@ -364,7 +348,7 @@ class P2PMatrixBuilder(MatrixBuilder): raise NotImplementedError("layer potentials on non-variables") kernel = expr.kernel.get_base_kernel() - kernel_args = _get_kernel_args(self, kernel, expr, source) + kernel_args = _get_kernel_args(self, kernel, expr, self.places.lpot) if self.exclude_self: kernel_args["target_to_source"] = \ cl.array.arange(self.queue, 0, target_discr.nnodes, dtype=np.int) @@ -388,17 +372,15 @@ class P2PMatrixBuilder(MatrixBuilder): # {{{ block matrix builders class MatrixBlockBuilderBase(EvaluationMapperBase): - def __init__(self, queue, dep_expr, other_dep_exprs, dep_source, + def __init__(self, queue, dep_expr, other_dep_exprs, dep_discr, places, context, index_set): super(MatrixBlockBuilderBase, self).__init__(context=context) self.queue = queue self.dep_expr = dep_expr self.other_dep_exprs = other_dep_exprs - self.dep_source = dep_source - self.dep_discr = dep_source.density_discr + self.dep_discr = dep_discr self.places = places - self.index_set = index_set def _map_dep_variable(self): @@ -448,7 +430,7 @@ class MatrixBlockBuilderBase(EvaluationMapperBase): rec_arg = cl.array.to_device(self.queue, rec_arg) op = expr.function(sym.var("u")) - result = bind(self.dep_source, op)(self.queue, u=rec_arg) + result = bind(self.place.lpot, op)(self.queue, u=rec_arg) if isinstance(result, cl.array.Array): result = result.get() @@ -457,13 +439,13 @@ class MatrixBlockBuilderBase(EvaluationMapperBase): class NearFieldBlockBuilder(MatrixBlockBuilderBase): - def __init__(self, queue, dep_expr, other_dep_exprs, dep_source, + def __init__(self, queue, dep_expr, other_dep_exprs, dep_discr, places, context, index_set): super(NearFieldBlockBuilder, self).__init__(queue, - dep_expr, other_dep_exprs, dep_source, places, context, index_set) + dep_expr, other_dep_exprs, dep_discr, places, context, index_set) self.dummy = MatrixBlockBuilderBase(queue, - dep_expr, other_dep_exprs, dep_source, places, context, index_set) + dep_expr, other_dep_exprs, dep_discr, places, context, index_set) def _map_dep_variable(self): tgtindices = self.index_set.row.indices.get(self.queue).reshape(-1, 1) @@ -472,17 +454,8 @@ class NearFieldBlockBuilder(MatrixBlockBuilderBase): return np.equal(tgtindices, srcindices).astype(np.float64) def map_int_g(self, expr): - where_source = expr.source - if where_source is sym.DEFAULT_SOURCE: - where_source = sym.QBXSourceStage1(where_source) - - where_target = expr.target - if where_target is sym.DEFAULT_TARGET: - where_target = sym.QBXSourceStage1(expr.target) - - from pytential.symbolic.execution import _get_discretization - source, source_discr = _get_discretization(self.places, where_source) - _, target_discr = _get_discretization(self.places, where_target) + source_discr = self.places[expr.source] + target_discr = self.places[expr.target] if source_discr is not target_discr: raise NotImplementedError() @@ -496,10 +469,10 @@ class NearFieldBlockBuilder(MatrixBlockBuilderBase): raise NotImplementedError("layer potentials on non-variables") kernel = expr.kernel - kernel_args = _get_layer_potential_args(self, expr, source) + kernel_args = _get_layer_potential_args(self, expr, self.places.lpot) from sumpy.expansion.local import LineTaylorLocalExpansion - local_expn = LineTaylorLocalExpansion(kernel, source.qbx_order) + local_expn = LineTaylorLocalExpansion(kernel, self.places.lpot.qbx_order) from sumpy.qbx import LayerPotentialMatrixBlockGenerator mat_gen = LayerPotentialMatrixBlockGenerator( @@ -507,7 +480,7 @@ class NearFieldBlockBuilder(MatrixBlockBuilderBase): assert abs(expr.qbx_forced_limit) > 0 centers, radii = _get_centers_and_expansion_radii(self.queue, - source, target_discr, expr.qbx_forced_limit) + self.places.lpot, target_discr, expr.qbx_forced_limit) _, (mat,) = mat_gen(self.queue, targets=target_discr.nodes(), @@ -517,7 +490,8 @@ class NearFieldBlockBuilder(MatrixBlockBuilderBase): index_set=self.index_set, **kernel_args) - waa = _get_weights_and_area_elements(self.queue, source, where_source) + waa = _get_weights_and_area_elements(self.queue, + self.places.lpot, source_discr) mat *= waa[self.index_set.linear_col_indices] mat = mat.get(self.queue) @@ -527,13 +501,13 @@ class NearFieldBlockBuilder(MatrixBlockBuilderBase): class FarFieldBlockBuilder(MatrixBlockBuilderBase): - def __init__(self, queue, dep_expr, other_dep_exprs, dep_source, + def __init__(self, queue, dep_expr, other_dep_exprs, dep_discr, places, context, index_set, exclude_self=True): super(FarFieldBlockBuilder, self).__init__(queue, - dep_expr, other_dep_exprs, dep_source, places, context, index_set) + dep_expr, other_dep_exprs, dep_discr, places, context, index_set) self.dummy = MatrixBlockBuilderBase(queue, - dep_expr, other_dep_exprs, dep_source, places, context, index_set) + dep_expr, other_dep_exprs, dep_discr, places, context, index_set) self.exclude_self = exclude_self def _map_dep_variable(self): @@ -543,20 +517,8 @@ class FarFieldBlockBuilder(MatrixBlockBuilderBase): return np.equal(tgtindices, srcindices).astype(np.float64) def map_int_g(self, expr): - where_source = expr.source - if where_source is sym.DEFAULT_SOURCE: - where_source = sym.QBXSourceStage1(where_source) - - where_target = expr.target - if where_target is sym.DEFAULT_TARGET: - where_target = sym.QBXSourceStage1(expr.target) - - from pytential.symbolic.execution import _get_discretization - source, source_discr = _get_discretization(self.places, where_source) - _, target_discr = _get_discretization(self.places, where_target) - - if source_discr is not target_discr: - raise NotImplementedError() + source_discr = self.places[expr.source] + target_discr = self.places[expr.target] rec_density = self.dummy.rec(expr.density) if is_zero(rec_density): @@ -567,7 +529,7 @@ class FarFieldBlockBuilder(MatrixBlockBuilderBase): raise NotImplementedError("layer potentials on non-variables") kernel = expr.kernel.get_base_kernel() - kernel_args = _get_kernel_args(self, kernel, expr, source) + kernel_args = _get_kernel_args(self, kernel, expr, self.places.lpot) if self.exclude_self: kernel_args["target_to_source"] = \ cl.array.arange(self.queue, 0, target_discr.nnodes, dtype=np.int) diff --git a/test/test_matrix.py b/test/test_matrix.py index 842029b1..58a7ea71 100644 --- a/test/test_matrix.py +++ b/test/test_matrix.py @@ -47,9 +47,23 @@ from pyopencl.tools import ( # noqa as pytest_generate_tests) -def _build_op(lpot_id, k=0, ndim=2): +def _build_op(lpot_id, + k=0, + ndim=2, + qbx_forced_limit=-1, + source=None, + target=None): + if source is DEFAULT_SOURCE: + source = QBXSourceQuadStage2(source) + if target is DEFAULT_TARGET: + target = QBXSourceStage1(target) + knl_kwargs = { + "qbx_forced_limit": qbx_forced_limit, + "source": source, + "target": target + } + from sumpy.kernel import LaplaceKernel, HelmholtzKernel - knl_kwargs = {"qbx_forced_limit": "avg"} if k: knl = HelmholtzKernel(ndim) knl_kwargs["k"] = k @@ -171,7 +185,7 @@ def test_matrix_build(ctx_factory, k, curve_f, lpot_id, visualize=False): abs_err = la.norm(res_mat - res_matvec, np.inf) rel_err = abs_err / la.norm(res_matvec, np.inf) - print(abs_err, rel_err) + print("AbsErr {:.5e} RelErr {:.5e}".format(abs_err, rel_err)) assert rel_err < 1e-13 @@ -197,9 +211,9 @@ def test_p2p_block_builder(ctx_factory, factor, ndim, lpot_id, tgtindices = _build_block_index(qbx.density_discr, method='nodes', factor=factor) - from pytential.symbolic.execution import prepare_places, prepare_expr - places = prepare_places(qbx) - expr = prepare_expr(places, op) + from pytential.symbolic.execution import BindingLocation, _prepare_expr + places = BindingLocation(qbx) + expr = _prepare_expr(places, op) from sumpy.tools import MatrixBlockIndexRanges index_set = MatrixBlockIndexRanges(ctx, tgtindices, srcindices) @@ -208,7 +222,7 @@ def test_p2p_block_builder(ctx_factory, factor, ndim, lpot_id, mbuilder = FarFieldBlockBuilder(queue, dep_expr=u_sym, other_dep_exprs=[], - dep_source=places[DEFAULT_SOURCE], + dep_discr=places[DEFAULT_SOURCE], places=places, context={}, index_set=index_set) @@ -218,8 +232,7 @@ def test_p2p_block_builder(ctx_factory, factor, ndim, lpot_id, mbuilder = P2PMatrixBuilder(queue, dep_expr=u_sym, other_dep_exprs=[], - dep_source=places[DEFAULT_SOURCE], - dep_discr=places[DEFAULT_SOURCE].density_discr, + dep_discr=places[DEFAULT_SOURCE], places=places, context={}) mat = mbuilder(expr) @@ -272,11 +285,14 @@ def test_qbx_block_builder(ctx_factory, ndim, lpot_id, visualize=False): tgtindices = _build_block_index(qbx.density_discr) srcindices = _build_block_index(qbx.density_discr) - from pytential.symbolic.execution import prepare_places, prepare_expr - where = (QBXSourceStage1(DEFAULT_SOURCE), - QBXSourceStage1(DEFAULT_TARGET)) - places = prepare_places(qbx, auto_where=where) - expr = prepare_expr(places, op, auto_where=where) + # NOTE: NearFieldBlockBuilder only does stage1/stage1 or stage2/stage2, + # so we need to hardcode the discr for MatrixBuilder too, since the + # defaults are different + where = (QBXSourceStage1(DEFAULT_SOURCE), QBXSourceStage1(DEFAULT_TARGET)) + + from pytential.symbolic.execution import BindingLocation, _prepare_expr + places = BindingLocation(qbx, auto_where=where) + expr = _prepare_expr(places, op) from sumpy.tools import MatrixBlockIndexRanges index_set = MatrixBlockIndexRanges(ctx, tgtindices, srcindices) @@ -285,19 +301,17 @@ def test_qbx_block_builder(ctx_factory, ndim, lpot_id, visualize=False): mbuilder = NearFieldBlockBuilder(queue, dep_expr=u_sym, other_dep_exprs=[], - dep_source=places[where[0]], + dep_discr=places[where[0]], places=places, context={}, index_set=index_set) blk = mbuilder(expr) - from pytential.symbolic.execution import _get_discretization from pytential.symbolic.matrix import MatrixBuilder mbuilder = MatrixBuilder(queue, dep_expr=u_sym, other_dep_exprs=[], - dep_source=places[where[0]], - dep_discr=_get_discretization(places, where[0])[1], + dep_discr=places[where[0]], places=places, context={}) mat = mbuilder(expr) @@ -349,8 +363,11 @@ def test_build_matrix_where(ctx_factory, where, visualize=False): from test_linalg_proxy import _build_qbx_discr qbx = _build_qbx_discr(queue, nelements=8, target_order=2, ndim=2, curve_f=partial(ellipse, 1.0)) - op, u_sym, _ = _build_op(lpot_id=1, ndim=2) + qbx_forced_limit = -1 + op, u_sym, _ = _build_op(lpot_id=1, ndim=2, + source=where[0], target=where[1], + qbx_forced_limit=qbx_forced_limit) # associate quad_stage2 targets to centers from pytential.qbx.target_assoc import associate_targets_to_qbx_centers @@ -393,10 +410,10 @@ def test_build_matrix_where(ctx_factory, where, visualize=False): from pytential.symbolic.execution import build_matrix mat = build_matrix(queue, qbx, op, u_sym, auto_where=where) - from pytential.symbolic.execution import _get_discretization - source_where, target_where = where - _, source_discr = _get_discretization(qbx, source_where) - _, target_discr = _get_discretization(qbx, target_where) + from pytential.symbolic.execution import BindingLocation + places = BindingLocation(qbx, auto_where=where) + source_discr = places[where[0]] + target_discr = places[where[1]] assert mat.shape == (target_discr.nnodes, source_discr.nnodes) -- GitLab From 9577be8faec5b3768eddade68ec99b704c8f2a79 Mon Sep 17 00:00:00 2001 From: Alex Fikl Date: Sun, 5 Aug 2018 20:07:39 -0500 Subject: [PATCH 19/34] flake8 and doc fixes --- pytential/symbolic/execution.py | 13 +++++-------- pytential/symbolic/matrix.py | 3 ++- 2 files changed, 7 insertions(+), 9 deletions(-) diff --git a/pytential/symbolic/execution.py b/pytential/symbolic/execution.py index 0217be3d..d6d2a88d 100644 --- a/pytential/symbolic/execution.py +++ b/pytential/symbolic/execution.py @@ -286,7 +286,7 @@ class MatVecOp: def _prepare_domains(nresults, places, domains, default_val): """ :arg nresults: number of results. - :arg places: result of :func:`prepare_places`. + :arg places: :class:`pytential.symbolic.execution.BindingLocation`. :arg domains: recommended domains. :arg default_val: default value for domains which are not provided. @@ -311,13 +311,10 @@ def _prepare_domains(nresults, places, domains, default_val): def _prepare_expr(places, expr): """ - :arg places: result of :func:`prepare_places`. - :arg expr: an array of symbolic expressions. - :arg auto_where: identifiers for source and/or target locations. If `None`, - `where` attributes are automatically found. - - :return: processed symbolic expressions, tagger with the given `where` - identifiers. + :arg places: :class:`pytential.symbolic.execution.BindingLocation`. + :arg expr: a symbolic expression. + :return: processed symbolic expressions, tagged with the appropriate + `where` identifier from places, etc. """ from pytential.source import LayerPotentialSourceBase diff --git a/pytential/symbolic/matrix.py b/pytential/symbolic/matrix.py index 8c9a354e..3e9a6c9d 100644 --- a/pytential/symbolic/matrix.py +++ b/pytential/symbolic/matrix.py @@ -271,7 +271,8 @@ class MatrixBuilder(EvaluationMapperBase): **kernel_args) mat = mat.get() - waa = _get_weights_and_area_elements(self.queue, self.places.lpot, source_discr) + waa = _get_weights_and_area_elements(self.queue, + self.places.lpot, source_discr) mat[:, :] *= waa.get(self.queue) if target_discr.nnodes != source_discr.nnodes: -- GitLab From ccdcb980b34b2ffdbed577518c5c0dea3b7434b5 Mon Sep 17 00:00:00 2001 From: Alex Fikl Date: Thu, 9 Aug 2018 20:26:52 -0500 Subject: [PATCH 20/34] execution: more work on making the geometry collection ready. * renamed to GeometryCollection * the `sources` are always a LayerPotentialSourceBase of some sort. * added a get_discretization that gets the actual mesh and used that in most places. * added process_optemplate to the PotentialSource base class, to match the documentation. --- pytential/source.py | 34 ++++-- pytential/symbolic/compiler.py | 4 +- pytential/symbolic/execution.py | 186 +++++++++++++++----------------- pytential/symbolic/mappers.py | 4 +- pytential/symbolic/matrix.py | 90 +++++++++------- test/test_matrix.py | 69 ++++++------ 6 files changed, 198 insertions(+), 189 deletions(-) diff --git a/pytential/source.py b/pytential/source.py index 6420494e..6110972e 100644 --- a/pytential/source.py +++ b/pytential/source.py @@ -38,7 +38,7 @@ __doc__ = """ class PotentialSource(object): """ - .. method:: preprocess_optemplate(name, expr) + .. method:: preprocess_optemplate(name, discretizations, expr) .. method:: op_group_features(expr) @@ -49,29 +49,43 @@ class PotentialSource(object): :class:`pytential.symbolic.primitives.IntG`. """ + def preprocess_optemplate(name, discretizations, expr): + return expr + # {{{ point potential source class PointPotentialSource(PotentialSource): """ - .. attribute:: points + .. attribute:: nodes - An :class:`pyopencl.array.Array` of shape ``[ambient_dim, npoints]``. + An :class:`pyopencl.array.Array` of shape ``[ambient_dim, nnodes]``. .. attribute:: nnodes """ - def __init__(self, cl_context, points): + def __init__(self, cl_context, nodes): self.cl_context = cl_context - self.points = points + self._nodes = nodes + + @property + def points(self): + from warnings import warn + warn("'points' has been renamed to nodes().", + DeprecationWarning, stacklevel=2) + + return self._nodes + + def nodes(self): + return self._nodes @property def real_dtype(self): - return self.points.dtype + return self._nodes.dtype @property def nnodes(self): - return self.points.shape[-1] + return self._nodes.shape[-1] @property def complex_dtype(self): @@ -82,7 +96,7 @@ class PointPotentialSource(PotentialSource): @property def ambient_dim(self): - return self.points.shape[0] + return self._nodes.shape[0] def op_group_features(self, expr): from sumpy.kernel import AxisTargetDerivativeRemover @@ -129,7 +143,7 @@ class PointPotentialSource(PotentialSource): p2p = self.get_p2p(insn.kernels) evt, output_for_each_kernel = p2p(queue, - target_discr.nodes(), self.points, + target_discr.nodes(), self._nodes, [strengths], **kernel_args) result.append((o.name, output_for_each_kernel[o.kernel_index])) @@ -139,7 +153,7 @@ class PointPotentialSource(PotentialSource): @memoize_method def weights_and_area_elements(self): with cl.CommandQueue(self.cl_context) as queue: - result = cl.array.empty(queue, self.points.shape[-1], + result = cl.array.empty(queue, self._nodes.shape[-1], dtype=self.real_dtype) result.fill(1) diff --git a/pytential/symbolic/compiler.py b/pytential/symbolic/compiler.py index ecb1fdca..17d23ebf 100644 --- a/pytential/symbolic/compiler.py +++ b/pytential/symbolic/compiler.py @@ -224,7 +224,7 @@ class ComputePotentialInstruction(Instruction): ", ".join(args), "\n ".join(lines)) def get_exec_function(self, exec_mapper): - source = exec_mapper.bound_expr.places.lpot + source = exec_mapper.bound_expr.places[self.source] return source.exec_compute_potential_insn def __hash__(self): @@ -444,7 +444,7 @@ class OperatorCompiler(IdentityMapper): def op_group_features(self, expr): from pytential.symbolic.primitives import hashable_kernel_args return ( - self.places.lpot.op_group_features(expr) + self.places[expr.source].op_group_features(expr) + hashable_kernel_args(expr.kernel_arguments)) @memoize_method diff --git a/pytential/symbolic/execution.py b/pytential/symbolic/execution.py index d6d2a88d..a83301fa 100644 --- a/pytential/symbolic/execution.py +++ b/pytential/symbolic/execution.py @@ -283,48 +283,49 @@ class MatVecOp: # {{{ expression prep -def _prepare_domains(nresults, places, domains, default_val): +def _prepare_domains(nresults, places, domains, default_domain): """ :arg nresults: number of results. - :arg places: :class:`pytential.symbolic.execution.BindingLocation`. + :arg places: :class:`pytential.symbolic.execution.GeometryCollection`. :arg domains: recommended domains. - :arg default_val: default value for domains which are not provided. + :arg default_domain: default value for domains which are not provided. :return: a list of domains for each result. If domains is `None`, each - element in the list is *default_val*. If *domains* is not a list - of domains, each element in the resulting list is *domains*. Otherwise, - *domains* is returned as is. + element in the list is *default_domain*. If *domains* is a scalar + (i.e., not a *list* or *tuple*), each element in the list is + *domains*. Otherwise, *domains* is returned as is. """ if domains is None: - if default_val not in places: + if default_domain not in places: raise RuntimeError("'domains is None' requires " - "default domain to be defined") - dom_name = default_val + "default domain to be defined in places") + dom_name = default_domain return nresults * [dom_name] - elif not isinstance(domains, (list, tuple)): + if not isinstance(domains, (list, tuple)): dom_name = domains return nresults * [dom_name] - else: - assert len(domains) == nresults - return domains + + assert len(domains) == nresults + return domains def _prepare_expr(places, expr): """ - :arg places: :class:`pytential.symbolic.execution.BindingLocation`. + :arg places: :class:`pytential.symbolic.execution.GeometryCollection`. :arg expr: a symbolic expression. :return: processed symbolic expressions, tagged with the appropriate `where` identifier from places, etc. """ - from pytential.source import LayerPotentialSourceBase + from pytential.source import PotentialSource from pytential.symbolic.mappers import ( ToTargetTagger, DerivativeBinder) expr = ToTargetTagger(*places.where)(expr) expr = DerivativeBinder()(expr) - if isinstance(places.lpot, LayerPotentialSourceBase): - expr = places.lpot.preprocess_optemplate(places.source, places, expr) + for name, place in six.iteritems(places.places): + if isinstance(places, PotentialSource): + expr = p.preprocess_optemplate(name, places, expr) return expr @@ -333,7 +334,7 @@ def _prepare_expr(places, expr): # {{{ bound expression -class BindingLocation(object): +class GeometryCollection(object): def __init__(self, places, auto_where=None): from meshmode.discretization import Discretization from pytential.source import LayerPotentialSourceBase, PotentialSource @@ -343,31 +344,24 @@ class BindingLocation(object): auto_where = DEFAULT_SOURCE, DEFAULT_TARGET self.where = auto_where - self.lpot = None self.places = {} if isinstance(places, LayerPotentialSourceBase): - self.lpot = places - self.places[self.source] = self[self.source] - self.places[self.target] = self[self.target] + self.places[self.source] = places + self.places[self.target] = \ + self.get_discretization(self.target, lpot=places) elif isinstance(places, (Discretization, TargetBase)): self.places[self.target] = places elif isinstance(places, tuple): source_discr, target_discr = places - if isinstance(source_discr, PotentialSource): - self.lpot = source_discr - if isinstance(source_discr, LayerPotentialSourceBase): - source_discr = self[self.source] - self.places[self.source] = source_discr self.places[self.target] = target_discr else: - for key, value in six.iteritems(places): - self.places[key] = value + self.places = places.copy() for p in six.itervalues(self.places): if not isinstance(p, (PotentialSource, TargetBase, Discretization)): raise TypeError("Must pass discretization, targets or " - "layer potential sources as 'places'") + "layer potential sources as 'places'.") self.caches = {} @@ -379,51 +373,38 @@ class BindingLocation(object): def target(self): return self.where[1] - def _get_lpot_discretization(self, where): - from pytential.source import LayerPotentialSourceBase - if not isinstance(self.lpot, LayerPotentialSourceBase): - return None - + def get_discretization(self, where, lpot=None): if where is DEFAULT_SOURCE or where is DEFAULT_TARGET: where = QBXSourceStage1(where) + if lpot is None: + if where in self.places: + lpot = self.places[where] + else: + lpot = self.places.get(where.where, None) + + from pytential.source import LayerPotentialSourceBase + if not isinstance(lpot, LayerPotentialSourceBase): + return lpot + if isinstance(where, QBXSourceStage1): - return self.lpot.density_discr + return lpot.density_discr if isinstance(where, QBXSourceStage2): - return self.lpot.stage2_density_discr + return lpot.stage2_density_discr if isinstance(where, QBXSourceQuadStage2): - return self.lpot.quad_stage2_density_discr + return lpot.quad_stage2_density_discr - return None + raise ValueError('Unknown `where` identifier: {}'.format(type(where))) def __getitem__(self, where): - if where in self.places: - discr = self.places[where] - else: - discr = self._get_lpot_discretization(where) - - if discr is None: - from pytential.symbolic.mappers import stringify_where - raise KeyError("Identifier '{}' not in places".format( - stringify_where(where))) - - return discr - - def __iter__(self): - return iter(self.places.keys()) - - def keys(self): - return self.places.keys() + return self.places[where] - def values(self): - return self.places.values() + def __contains__(self, where): + return where in self.places def items(self): return self.places.items() - def __hash__(self): - return hash(tuple(self.values())) - def get_cache(self, name): return self.caches.setdefault(name, {}) @@ -441,7 +422,7 @@ class BoundExpression(object): return self.places.get_cache(name) def get_discretization(self, where): - return self.places[where] + return self.places.get_discretization(where) def scipy_op(self, queue, arg_name, dtype, domains=None, **extra_args): """ @@ -488,14 +469,15 @@ class BoundExpression(object): def bind(places, expr, auto_where=None): """ - :arg places: a mapping of symbolic names to - :class:`pytential.discretization.Discretization` objects or a subclass - of :class:`pytential.discretization.target.TargetBase`. - :arg auto_where: For simple source-to-self or source-to-target + :arg places: a collection or mapping of symbolic names to + :class:`meshmode.discretization.Discretization` objects, subclasses + of :class:`pytential.source.PotentialSource` or subclasses of + :class:`pytential.target.TargetBase`. + :arg auto_where: for simple source-to-self or source-to-target evaluations, find 'where' attributes automatically. """ - places = BindingLocation(places, auto_where=auto_where) + places = GeometryCollection(places, auto_where=auto_where) expr = _prepare_expr(places, expr) return BoundExpression(places, expr) @@ -505,14 +487,44 @@ def bind(places, expr, auto_where=None): # {{{ matrix building +def _bmat(blocks, dtypes): + from pytools import single_valued + from pytential.symbolic.matrix import is_zero + + nrows = blocks.shape[0] + ncolumns = blocks.shape[1] + + # "block row starts"/"block column starts" + brs = np.cumsum([0] + + [single_valued(blocks[ibrow, ibcol].shape[0] + for ibcol in range(ncolumns) + if not is_zero(blocks[ibrow, ibcol])) + for ibrow in range(nrows)]) + + bcs = np.cumsum([0] + + [single_valued(blocks[ibrow, ibcol].shape[1] + for ibrow in range(nrows) + if not is_zero(blocks[ibrow, ibcol])) + for ibcol in range(ncolumns)]) + + result = np.zeros((brs[-1], bcs[-1]), + dtype=np.find_common_type(dtypes, [])) + for ibcol in range(ncolumns): + for ibrow in range(nrows): + result[brs[ibrow]:brs[ibrow + 1], bcs[ibcol]:bcs[ibcol + 1]] = \ + blocks[ibrow, ibcol] + + return result + + def build_matrix(queue, places, exprs, input_exprs, domains=None, auto_where=None, context=None): """ :arg queue: a :class:`pyopencl.CommandQueue` used to synchronize the calculation. - :arg places: a mapping of symbolic names to - :class:`pytential.discretization.Discretization` objects or a subclass - of :class:`pytential.discretization.target.TargetBase`. + :arg places: a collection or mapping of symbolic names to + :class:`meshmode.discretization.Discretization` objects, subclasses + of :class:`pytential.source.PotentialSource` or subclasses of :arg input_exprs: An object array of expressions corresponding to the input block columns of the matrix. @@ -521,7 +533,7 @@ def build_matrix(queue, places, exprs, input_exprs, domains=None, *None* values indicating the domains on which each component of the solution vector lives. *None* values indicate that the component is a scalar. If *None*, - :class:`pytential.symbolic.primitives.DEFAULT_TARGET`, is required + :class:`pytential.symbolic.primitives.DEFAULT_SOURCE`, is required to be a key in :attr:`places`. :arg auto_where: For simple source-to-self or source-to-target evaluations, find 'where' attributes automatically. @@ -531,7 +543,7 @@ def build_matrix(queue, places, exprs, input_exprs, domains=None, context = {} from pytools.obj_array import is_obj_array, make_obj_array - places = BindingLocation(places, auto_where=auto_where) + places = GeometryCollection(places, auto_where=auto_where) exprs = _prepare_expr(places, exprs) if not is_obj_array(exprs): @@ -556,7 +568,8 @@ def build_matrix(queue, places, exprs, input_exprs, domains=None, dep_expr=input_exprs[ibcol], other_dep_exprs=(input_exprs[:ibcol] + input_exprs[ibcol + 1:]), - dep_discr=places[domains[ibcol]], + dep_source=places[domains[ibcol]], + dep_discr=places.get_discretization(domains[ibcol]), places=places, context=context) @@ -568,32 +581,7 @@ def build_matrix(queue, places, exprs, input_exprs, domains=None, if isinstance(block, np.ndarray): dtypes.append(block.dtype) - from pytools import single_valued - block_row_counts = [ - single_valued( - blocks[ibrow, ibcol].shape[0] - for ibcol in range(nblock_columns) - if not is_zero(blocks[ibrow, ibcol])) - for ibrow in range(nblock_rows)] - - block_col_counts = [ - single_valued( - blocks[ibrow, ibcol].shape[1] - for ibrow in range(nblock_rows) - if not is_zero(blocks[ibrow, ibcol])) - for ibcol in range(nblock_columns)] - - # "block row starts"/"block column starts" - brs = np.cumsum([0] + block_row_counts) - bcs = np.cumsum([0] + block_col_counts) - - result = np.zeros((brs[-1], bcs[-1]), dtype=np.find_common_type(dtypes, [])) - for ibcol in range(nblock_columns): - for ibrow in range(nblock_rows): - result[brs[ibrow]:brs[ibrow + 1], bcs[ibcol]:bcs[ibcol + 1]] = \ - blocks[ibrow, ibcol] - - return cl.array.to_device(queue, result) + return cl.array.to_device(queue, _bmat(blocks, dtypes)) # }}} diff --git a/pytential/symbolic/mappers.py b/pytential/symbolic/mappers.py index a8696efe..00c43689 100644 --- a/pytential/symbolic/mappers.py +++ b/pytential/symbolic/mappers.py @@ -399,8 +399,8 @@ class QBXPreprocessor(IdentityMapper): self.places = places def map_int_g(self, expr): - source_discr = self.places[self.source_name] - target_discr = self.places[expr.target] + source_discr = self.places.get_discretization(self.source_name) + target_discr = self.places.get_discretization(expr.target) if expr.qbx_forced_limit == 0: raise ValueError("qbx_forced_limit == 0 was a bad idea and " diff --git a/pytential/symbolic/matrix.py b/pytential/symbolic/matrix.py index 3e9a6c9d..79db0f8d 100644 --- a/pytential/symbolic/matrix.py +++ b/pytential/symbolic/matrix.py @@ -41,7 +41,8 @@ def is_zero(x): def _resample_arg(queue, source, x): - if source is None: + from pytential.source import LayerPotentialSourceBase + if not isinstance(source, LayerPotentialSourceBase): return x if not isinstance(x, np.ndarray): @@ -138,13 +139,14 @@ def _get_centers_and_expansion_radii(queue, source, target_discr, qbx_forced_lim # We'll cheat and build the matrix on the host. class MatrixBuilder(EvaluationMapperBase): - def __init__(self, queue, dep_expr, other_dep_exprs, dep_discr, - places, context): + def __init__(self, queue, dep_expr, other_dep_exprs, + dep_source, dep_discr, places, context): super(MatrixBuilder, self).__init__(context=context) self.queue = queue self.dep_expr = dep_expr self.other_dep_exprs = other_dep_exprs + self.dep_source = dep_source self.dep_discr = dep_discr self.places = places @@ -233,13 +235,14 @@ class MatrixBuilder(EvaluationMapperBase): return vecs_and_scalars def map_int_g(self, expr): - # TODO: should this go into QBXPreprocessor / LocationTagger + # TODO: should this go into QBXPreprocessor / LocationTagger? where_source = expr.source if where_source is sym.DEFAULT_SOURCE: where_source = sym.QBXSourceQuadStage2(expr.source) - source_discr = self.places[where_source] - target_discr = self.places[expr.target] + source = self.places[expr.source] + source_discr = self.places.get_discretization(where_source) + target_discr = self.places.get_discretization(expr.target) rec_density = self.rec(expr.density) if is_zero(rec_density): @@ -250,10 +253,10 @@ class MatrixBuilder(EvaluationMapperBase): raise NotImplementedError("layer potentials on non-variables") kernel = expr.kernel - kernel_args = _get_layer_potential_args(self, expr, self.places.lpot) + kernel_args = _get_layer_potential_args(self, expr, source) from sumpy.expansion.local import LineTaylorLocalExpansion - local_expn = LineTaylorLocalExpansion(kernel, self.places.lpot.qbx_order) + local_expn = LineTaylorLocalExpansion(kernel, source.qbx_order) from sumpy.qbx import LayerPotentialMatrixGenerator mat_gen = LayerPotentialMatrixGenerator( @@ -261,7 +264,7 @@ class MatrixBuilder(EvaluationMapperBase): assert abs(expr.qbx_forced_limit) > 0 centers, radii = _get_centers_and_expansion_radii(self.queue, - self.places.lpot, target_discr, expr.qbx_forced_limit) + source, target_discr, expr.qbx_forced_limit) _, (mat,) = mat_gen(self.queue, targets=target_discr.nodes(), @@ -271,14 +274,13 @@ class MatrixBuilder(EvaluationMapperBase): **kernel_args) mat = mat.get() - waa = _get_weights_and_area_elements(self.queue, - self.places.lpot, source_discr) + waa = _get_weights_and_area_elements(self.queue, source, source_discr) mat[:, :] *= waa.get(self.queue) if target_discr.nnodes != source_discr.nnodes: assert target_discr.nnodes < source_discr.nnodes - resampler = self.places.lpot.direct_resampler + resampler = source.direct_resampler resample_mat = resampler.full_resample_matrix(self.queue).get(self.queue) mat = mat.dot(resample_mat) @@ -316,7 +318,7 @@ class MatrixBuilder(EvaluationMapperBase): rec_arg = cl.array.to_device(self.queue, rec_arg) op = expr.function(sym.var("u")) - result = bind(self.places.lpot, op)(self.queue, u=rec_arg) + result = bind(self.dep_source, op)(self.queue, u=rec_arg) if isinstance(result, cl.array.Array): result = result.get() @@ -329,16 +331,18 @@ class MatrixBuilder(EvaluationMapperBase): # {{{ p2p matrix builder class P2PMatrixBuilder(MatrixBuilder): - def __init__(self, queue, dep_expr, other_dep_exprs, dep_discr, - places, context, exclude_self=True): + def __init__(self, queue, dep_expr, other_dep_exprs, + dep_source, dep_discr, places, context, exclude_self=True): super(P2PMatrixBuilder, self).__init__(queue, - dep_expr, other_dep_exprs, dep_discr, places, context) + dep_expr, other_dep_exprs, dep_source, dep_discr, + places, context) self.exclude_self = exclude_self def map_int_g(self, expr): - source_discr = self.places[expr.source] - target_discr = self.places[expr.target] + source = self.places[expr.source] + source_discr = self.places.get_discretization(expr.source) + target_discr = self.places.get_discretization(expr.target) rec_density = self.rec(expr.density) if is_zero(rec_density): @@ -349,7 +353,7 @@ class P2PMatrixBuilder(MatrixBuilder): raise NotImplementedError("layer potentials on non-variables") kernel = expr.kernel.get_base_kernel() - kernel_args = _get_kernel_args(self, kernel, expr, self.places.lpot) + kernel_args = _get_kernel_args(self, kernel, expr, source) if self.exclude_self: kernel_args["target_to_source"] = \ cl.array.arange(self.queue, 0, target_discr.nnodes, dtype=np.int) @@ -373,13 +377,14 @@ class P2PMatrixBuilder(MatrixBuilder): # {{{ block matrix builders class MatrixBlockBuilderBase(EvaluationMapperBase): - def __init__(self, queue, dep_expr, other_dep_exprs, dep_discr, - places, context, index_set): + def __init__(self, queue, dep_expr, other_dep_exprs, + dep_source, dep_discr, places, index_set, context): super(MatrixBlockBuilderBase, self).__init__(context=context) self.queue = queue self.dep_expr = dep_expr self.other_dep_exprs = other_dep_exprs + self.dep_source = dep_source self.dep_discr = dep_discr self.places = places self.index_set = index_set @@ -431,7 +436,7 @@ class MatrixBlockBuilderBase(EvaluationMapperBase): rec_arg = cl.array.to_device(self.queue, rec_arg) op = expr.function(sym.var("u")) - result = bind(self.place.lpot, op)(self.queue, u=rec_arg) + result = bind(self.dep_source, op)(self.queue, u=rec_arg) if isinstance(result, cl.array.Array): result = result.get() @@ -440,13 +445,15 @@ class MatrixBlockBuilderBase(EvaluationMapperBase): class NearFieldBlockBuilder(MatrixBlockBuilderBase): - def __init__(self, queue, dep_expr, other_dep_exprs, dep_discr, - places, context, index_set): + def __init__(self, queue, dep_expr, other_dep_exprs, dep_source, dep_discr, + places, index_set, context): super(NearFieldBlockBuilder, self).__init__(queue, - dep_expr, other_dep_exprs, dep_discr, places, context, index_set) + dep_expr, other_dep_exprs, dep_source, dep_discr, + places, index_set, context) self.dummy = MatrixBlockBuilderBase(queue, - dep_expr, other_dep_exprs, dep_discr, places, context, index_set) + dep_expr, other_dep_exprs, dep_source, dep_discr, + places, index_set, context) def _map_dep_variable(self): tgtindices = self.index_set.row.indices.get(self.queue).reshape(-1, 1) @@ -455,8 +462,9 @@ class NearFieldBlockBuilder(MatrixBlockBuilderBase): return np.equal(tgtindices, srcindices).astype(np.float64) def map_int_g(self, expr): - source_discr = self.places[expr.source] - target_discr = self.places[expr.target] + source = self.places[expr.source] + source_discr = self.places.get_discretization(expr.source) + target_discr = self.places.get_discretization(expr.target) if source_discr is not target_discr: raise NotImplementedError() @@ -470,10 +478,10 @@ class NearFieldBlockBuilder(MatrixBlockBuilderBase): raise NotImplementedError("layer potentials on non-variables") kernel = expr.kernel - kernel_args = _get_layer_potential_args(self, expr, self.places.lpot) + kernel_args = _get_layer_potential_args(self, expr, source) from sumpy.expansion.local import LineTaylorLocalExpansion - local_expn = LineTaylorLocalExpansion(kernel, self.places.lpot.qbx_order) + local_expn = LineTaylorLocalExpansion(kernel, source.qbx_order) from sumpy.qbx import LayerPotentialMatrixBlockGenerator mat_gen = LayerPotentialMatrixBlockGenerator( @@ -481,7 +489,7 @@ class NearFieldBlockBuilder(MatrixBlockBuilderBase): assert abs(expr.qbx_forced_limit) > 0 centers, radii = _get_centers_and_expansion_radii(self.queue, - self.places.lpot, target_discr, expr.qbx_forced_limit) + source, target_discr, expr.qbx_forced_limit) _, (mat,) = mat_gen(self.queue, targets=target_discr.nodes(), @@ -491,8 +499,7 @@ class NearFieldBlockBuilder(MatrixBlockBuilderBase): index_set=self.index_set, **kernel_args) - waa = _get_weights_and_area_elements(self.queue, - self.places.lpot, source_discr) + waa = _get_weights_and_area_elements(self.queue, source, source_discr) mat *= waa[self.index_set.linear_col_indices] mat = mat.get(self.queue) @@ -502,13 +509,15 @@ class NearFieldBlockBuilder(MatrixBlockBuilderBase): class FarFieldBlockBuilder(MatrixBlockBuilderBase): - def __init__(self, queue, dep_expr, other_dep_exprs, dep_discr, - places, context, index_set, exclude_self=True): + def __init__(self, queue, dep_expr, other_dep_exprs, dep_source, dep_discr, + places, index_set, context, exclude_self=True): super(FarFieldBlockBuilder, self).__init__(queue, - dep_expr, other_dep_exprs, dep_discr, places, context, index_set) + dep_expr, other_dep_exprs, dep_source, dep_discr, + places, index_set, context) self.dummy = MatrixBlockBuilderBase(queue, - dep_expr, other_dep_exprs, dep_discr, places, context, index_set) + dep_expr, other_dep_exprs, dep_source, dep_discr, + places, index_set, context) self.exclude_self = exclude_self def _map_dep_variable(self): @@ -518,8 +527,9 @@ class FarFieldBlockBuilder(MatrixBlockBuilderBase): return np.equal(tgtindices, srcindices).astype(np.float64) def map_int_g(self, expr): - source_discr = self.places[expr.source] - target_discr = self.places[expr.target] + source = self.places[expr.source] + source_discr = self.places.get_discretization(expr.source) + target_discr = self.places.get_discretization(expr.target) rec_density = self.dummy.rec(expr.density) if is_zero(rec_density): @@ -530,7 +540,7 @@ class FarFieldBlockBuilder(MatrixBlockBuilderBase): raise NotImplementedError("layer potentials on non-variables") kernel = expr.kernel.get_base_kernel() - kernel_args = _get_kernel_args(self, kernel, expr, self.places.lpot) + kernel_args = _get_kernel_args(self, kernel, expr, source) if self.exclude_self: kernel_args["target_to_source"] = \ cl.array.arange(self.queue, 0, target_discr.nnodes, dtype=np.int) diff --git a/test/test_matrix.py b/test/test_matrix.py index 58a7ea71..487fcf8a 100644 --- a/test/test_matrix.py +++ b/test/test_matrix.py @@ -50,45 +50,37 @@ from pyopencl.tools import ( # noqa def _build_op(lpot_id, k=0, ndim=2, - qbx_forced_limit=-1, - source=None, - target=None): - if source is DEFAULT_SOURCE: - source = QBXSourceQuadStage2(source) - if target is DEFAULT_TARGET: - target = QBXSourceStage1(target) - knl_kwargs = { - "qbx_forced_limit": qbx_forced_limit, - "source": source, - "target": target - } + qbx_forced_limit=-1): from sumpy.kernel import LaplaceKernel, HelmholtzKernel if k: knl = HelmholtzKernel(ndim) - knl_kwargs["k"] = k + knl_kwargs = {"k": k} else: knl = LaplaceKernel(ndim) + knl_kwargs = {} + lpot_kwargs = {"qbx_forced_limit": qbx_forced_limit} + lpot_kwargs.update(knl_kwargs) if lpot_id == 1: # scalar single-layer potential u_sym = sym.var("u") - op = sym.S(knl, u_sym, **knl_kwargs) + op = sym.S(knl, u_sym, **lpot_kwargs) elif lpot_id == 2: # scalar double-layer potential u_sym = sym.var("u") - op = sym.D(knl, u_sym, **knl_kwargs) + op = sym.D(knl, u_sym, **lpot_kwargs) elif lpot_id == 3: # vector potential u_sym = sym.make_sym_vector("u", 2) u0_sym, u1_sym = u_sym op = make_obj_array([ - sym.Sp(knl, u0_sym, **knl_kwargs) + - sym.D(knl, u1_sym, **knl_kwargs), + sym.Sp(knl, u0_sym, **lpot_kwargs) + + sym.D(knl, u1_sym, **lpot_kwargs), - sym.S(knl, 0.4 * u0_sym, **knl_kwargs) + - 0.3 * sym.D(knl, u0_sym, **knl_kwargs) + sym.S(knl, 0.4 * u0_sym, **lpot_kwargs) + + 0.3 * sym.D(knl, u0_sym, **lpot_kwargs) ]) else: raise ValueError("Unknown lpot_id: {}".format(lpot_id)) @@ -211,9 +203,11 @@ def test_p2p_block_builder(ctx_factory, factor, ndim, lpot_id, tgtindices = _build_block_index(qbx.density_discr, method='nodes', factor=factor) - from pytential.symbolic.execution import BindingLocation, _prepare_expr - places = BindingLocation(qbx) + from pytential.symbolic.execution import GeometryCollection + from pytential.symbolic.execution import _prepare_expr, _prepare_domains + places = GeometryCollection(qbx) expr = _prepare_expr(places, op) + domains = _prepare_domains(1, places, None, places.source) from sumpy.tools import MatrixBlockIndexRanges index_set = MatrixBlockIndexRanges(ctx, tgtindices, srcindices) @@ -222,17 +216,19 @@ def test_p2p_block_builder(ctx_factory, factor, ndim, lpot_id, mbuilder = FarFieldBlockBuilder(queue, dep_expr=u_sym, other_dep_exprs=[], - dep_discr=places[DEFAULT_SOURCE], + dep_source=places[domains[0]], + dep_discr=places.get_discretization(domains[0]), places=places, - context={}, - index_set=index_set) + index_set=index_set, + context={}) blk = mbuilder(expr) from pytential.symbolic.matrix import P2PMatrixBuilder mbuilder = P2PMatrixBuilder(queue, dep_expr=u_sym, other_dep_exprs=[], - dep_discr=places[DEFAULT_SOURCE], + dep_source=places[domains[0]], + dep_discr=places.get_discretization(domains[0]), places=places, context={}) mat = mbuilder(expr) @@ -290,8 +286,8 @@ def test_qbx_block_builder(ctx_factory, ndim, lpot_id, visualize=False): # defaults are different where = (QBXSourceStage1(DEFAULT_SOURCE), QBXSourceStage1(DEFAULT_TARGET)) - from pytential.symbolic.execution import BindingLocation, _prepare_expr - places = BindingLocation(qbx, auto_where=where) + from pytential.symbolic.execution import GeometryCollection, _prepare_expr + places = GeometryCollection(qbx, auto_where=where) expr = _prepare_expr(places, op) from sumpy.tools import MatrixBlockIndexRanges @@ -301,17 +297,19 @@ def test_qbx_block_builder(ctx_factory, ndim, lpot_id, visualize=False): mbuilder = NearFieldBlockBuilder(queue, dep_expr=u_sym, other_dep_exprs=[], - dep_discr=places[where[0]], + dep_source=places[where[0]], + dep_discr=places.get_discretization(where[0]), places=places, - context={}, - index_set=index_set) + index_set=index_set, + context={}) blk = mbuilder(expr) from pytential.symbolic.matrix import MatrixBuilder mbuilder = MatrixBuilder(queue, dep_expr=u_sym, other_dep_exprs=[], - dep_discr=places[where[0]], + dep_source=places[where[0]], + dep_discr=places.get_discretization(where[0]), places=places, context={}) mat = mbuilder(expr) @@ -366,7 +364,6 @@ def test_build_matrix_where(ctx_factory, where, visualize=False): qbx_forced_limit = -1 op, u_sym, _ = _build_op(lpot_id=1, ndim=2, - source=where[0], target=where[1], qbx_forced_limit=qbx_forced_limit) # associate quad_stage2 targets to centers @@ -410,10 +407,10 @@ def test_build_matrix_where(ctx_factory, where, visualize=False): from pytential.symbolic.execution import build_matrix mat = build_matrix(queue, qbx, op, u_sym, auto_where=where) - from pytential.symbolic.execution import BindingLocation - places = BindingLocation(qbx, auto_where=where) - source_discr = places[where[0]] - target_discr = places[where[1]] + from pytential.symbolic.execution import GeometryCollection + places = GeometryCollection(qbx, auto_where=where) + source_discr = places.get_discretization(where[0]) + target_discr = places.get_discretization(where[1]) assert mat.shape == (target_discr.nnodes, source_discr.nnodes) -- GitLab From 5188a405013c7ff29deaa2b7be2d702aff457f02 Mon Sep 17 00:00:00 2001 From: Alex Fikl Date: Thu, 9 Aug 2018 20:57:54 -0500 Subject: [PATCH 21/34] execution: fix typo in prepare_expr --- pytential/source.py | 2 +- pytential/symbolic/execution.py | 7 ++++--- 2 files changed, 5 insertions(+), 4 deletions(-) diff --git a/pytential/source.py b/pytential/source.py index 6110972e..3bd58e46 100644 --- a/pytential/source.py +++ b/pytential/source.py @@ -49,7 +49,7 @@ class PotentialSource(object): :class:`pytential.symbolic.primitives.IntG`. """ - def preprocess_optemplate(name, discretizations, expr): + def preprocess_optemplate(self, name, discretizations, expr): return expr diff --git a/pytential/symbolic/execution.py b/pytential/symbolic/execution.py index a83301fa..db00c8ed 100644 --- a/pytential/symbolic/execution.py +++ b/pytential/symbolic/execution.py @@ -317,15 +317,16 @@ def _prepare_expr(places, expr): `where` identifier from places, etc. """ - from pytential.source import PotentialSource + from pytential.source import LayerPotentialSourceBase from pytential.symbolic.mappers import ( ToTargetTagger, DerivativeBinder) expr = ToTargetTagger(*places.where)(expr) expr = DerivativeBinder()(expr) + for name, place in six.iteritems(places.places): - if isinstance(places, PotentialSource): - expr = p.preprocess_optemplate(name, places, expr) + if isinstance(place, LayerPotentialSourceBase): + expr = place.preprocess_optemplate(name, places, expr) return expr -- GitLab From 451dffce306e9eb6afc65794350ce768ef9a387d Mon Sep 17 00:00:00 2001 From: Alex Fikl Date: Fri, 10 Aug 2018 10:16:59 -0500 Subject: [PATCH 22/34] matrix: fix a target_association bug and add more tests --- pytential/symbolic/execution.py | 10 ++-- pytential/symbolic/matrix.py | 6 +-- test/test_matrix.py | 85 ++++++++++++++++++--------------- 3 files changed, 55 insertions(+), 46 deletions(-) diff --git a/pytential/symbolic/execution.py b/pytential/symbolic/execution.py index db00c8ed..3eca37df 100644 --- a/pytential/symbolic/execution.py +++ b/pytential/symbolic/execution.py @@ -403,8 +403,8 @@ class GeometryCollection(object): def __contains__(self, where): return where in self.places - def items(self): - return self.places.items() + def copy(self): + return GeometryCollection(self.places, auto_where=self.where) def get_cache(self, name): return self.caches.setdefault(name, {}) @@ -478,7 +478,8 @@ def bind(places, expr, auto_where=None): evaluations, find 'where' attributes automatically. """ - places = GeometryCollection(places, auto_where=auto_where) + if not isinstance(places, GeometryCollection): + places = GeometryCollection(places, auto_where=auto_where) expr = _prepare_expr(places, expr) return BoundExpression(places, expr) @@ -544,7 +545,8 @@ def build_matrix(queue, places, exprs, input_exprs, domains=None, context = {} from pytools.obj_array import is_obj_array, make_obj_array - places = GeometryCollection(places, auto_where=auto_where) + if not isinstance(places, GeometryCollection): + places = GeometryCollection(places, auto_where=auto_where) exprs = _prepare_expr(places, exprs) if not is_obj_array(exprs): diff --git a/pytential/symbolic/matrix.py b/pytential/symbolic/matrix.py index 79db0f8d..f8c659c7 100644 --- a/pytential/symbolic/matrix.py +++ b/pytential/symbolic/matrix.py @@ -106,7 +106,7 @@ def _get_centers_and_expansion_radii(queue, source, target_discr, qbx_forced_lim else: from pytential.qbx.utils import get_interleaved_centers centers = get_interleaved_centers(queue, source) - radii = source._expansion_radii('nsources') + radii = source._expansion_radii('ncenters') # NOTE: using a very small tolerance to make sure all the stage2 # targets are associated to a center. We can't use the user provided @@ -124,9 +124,7 @@ def _get_centers_and_expansion_radii(queue, source, target_discr, qbx_forced_lim centers = [cl.array.take(c, assoc.target_to_center, queue=queue) for c in centers] - radii = cl.array.take(radii, - (assoc.target_to_center.with_queue(queue) / 2.0).astype(np.int), - queue=queue) + radii = cl.array.take(radii, assoc.target_to_center, queue=queue) return centers, radii diff --git a/test/test_matrix.py b/test/test_matrix.py index 487fcf8a..d16dd32a 100644 --- a/test/test_matrix.py +++ b/test/test_matrix.py @@ -366,44 +366,7 @@ def test_build_matrix_where(ctx_factory, where, visualize=False): op, u_sym, _ = _build_op(lpot_id=1, ndim=2, qbx_forced_limit=qbx_forced_limit) - # associate quad_stage2 targets to centers - from pytential.qbx.target_assoc import associate_targets_to_qbx_centers - - code_container = qbx.target_association_code_container - target_discr = qbx.quad_stage2_density_discr - target_assoc = associate_targets_to_qbx_centers( - qbx, - code_container.get_wrangler(queue), - [(target_discr, qbx_forced_limit)], - target_association_tolerance=1.0e-1).get(queue) - target_assoc = target_assoc.target_to_center - - if qbx.ambient_dim == 2 and visualize: - import matplotlib.pyplot as pt - from pytential.qbx.utils import get_interleaved_centers - sources = qbx.density_discr.nodes().get(queue) - targets = target_discr.nodes().get(queue) - centers = get_interleaved_centers(queue, qbx) - centers = np.vstack(c.get(queue) for c in centers) - radii = qbx._expansion_radii('nsources').get(queue) - - for i in range(centers[0].size): - itgt = np.where(target_assoc == i)[0] - if not len(itgt): - continue - - pt.figure(figsize=(10, 8), dpi=300) - pt.plot(sources[0], sources[1], 'k') - pt.plot(targets[0], targets[1], 'ko') - - line = pt.plot(targets[0, itgt], targets[1, itgt], 'o') - c = pt.Circle([centers[0][i], centers[1][i]], radii[i // 2], - color=line[0].get_color(), alpha=0.5) - pt.gca().add_artist(c) - - pt.savefig('test_assoc_quad_targets_to_centers_{:05}.png'.format(i // 2)) - pt.close() - + # build full QBX matrix from pytential.symbolic.execution import build_matrix mat = build_matrix(queue, qbx, op, u_sym, auto_where=where) @@ -414,6 +377,52 @@ def test_build_matrix_where(ctx_factory, where, visualize=False): assert mat.shape == (target_discr.nnodes, source_discr.nnodes) + # build full p2p matrix + from pytential.symbolic.execution import _prepare_expr + op = _prepare_expr(places, op) + + from pytential.symbolic.matrix import P2PMatrixBuilder + mbuilder = P2PMatrixBuilder(queue, + dep_expr=u_sym, + other_dep_exprs=[], + dep_source=places[where[0]], + dep_discr=places.get_discretization(where[0]), + places=places, + context={}) + mat = mbuilder(op) + + assert mat.shape == (target_discr.nnodes, source_discr.nnodes) + + # build block qbx and p2p matrices + from test_linalg_proxy import _build_block_index + srcindices = _build_block_index(source_discr, method='nodes', factor=0.6) + tgtindices = _build_block_index(target_discr, method='nodes', factor=0.6) + + from sumpy.tools import MatrixBlockIndexRanges + index_set = MatrixBlockIndexRanges(ctx, tgtindices, srcindices) + + from pytential.symbolic.matrix import NearFieldBlockBuilder + mbuilder = NearFieldBlockBuilder(queue, + dep_expr=u_sym, + other_dep_exprs=[], + dep_source=places[where[0]], + dep_discr=places.get_discretization(where[0]), + places=places, + index_set=index_set, + context={}) + mat = mbuilder(op) + + from pytential.symbolic.matrix import FarFieldBlockBuilder + mbuilder = FarFieldBlockBuilder(queue, + dep_expr=u_sym, + other_dep_exprs=[], + dep_source=places[where[0]], + dep_discr=places.get_discretization(where[0]), + places=places, + index_set=index_set, + context={}) + mat = mbuilder(op) + if __name__ == "__main__": import sys -- GitLab From defcbf55f0d473607bee31ab1fddbdf17b854172 Mon Sep 17 00:00:00 2001 From: Alex Fikl Date: Fri, 10 Aug 2018 10:39:26 -0500 Subject: [PATCH 23/34] matrix: add some docs --- pytential/symbolic/execution.py | 16 ++++++------ pytential/symbolic/matrix.py | 44 +++++++++++++++++++++++++++++++++ 2 files changed, 53 insertions(+), 7 deletions(-) diff --git a/pytential/symbolic/execution.py b/pytential/symbolic/execution.py index 3eca37df..b1e8eb4a 100644 --- a/pytential/symbolic/execution.py +++ b/pytential/symbolic/execution.py @@ -295,6 +295,7 @@ def _prepare_domains(nresults, places, domains, default_domain): (i.e., not a *list* or *tuple*), each element in the list is *domains*. Otherwise, *domains* is returned as is. """ + if domains is None: if default_domain not in places: raise RuntimeError("'domains is None' requires " @@ -522,15 +523,16 @@ def _bmat(blocks, dtypes): def build_matrix(queue, places, exprs, input_exprs, domains=None, auto_where=None, context=None): """ - :arg queue: a :class:`pyopencl.CommandQueue` used to synchronize - the calculation. - :arg places: a collection or mapping of symbolic names to + :arg queue: a :class:`pyopencl.CommandQueue`. + :arg places: a collection or mapping of symbolic names (see also + :class:`pytential.symbolic.execution.GeometryCollection`) to :class:`meshmode.discretization.Discretization` objects, subclasses of :class:`pytential.source.PotentialSource` or subclasses of - :arg input_exprs: An object array of expressions corresponding to the - input block columns of the matrix. - - May also be a single expression. + :class:`pytential.target.TargetBase`. + :arg exprs: an array of expressions corresponding to the output block + rows of the matrix. May also be a single expression. + :arg input_exprs: an array of expressions corresponding to the + input block columns of the matrix. May also be a single expression. :arg domains: a list of discretization identifiers (see 'places') or *None* values indicating the domains on which each component of the solution vector lives. *None* values indicate that the component diff --git a/pytential/symbolic/matrix.py b/pytential/symbolic/matrix.py index f8c659c7..16396f6c 100644 --- a/pytential/symbolic/matrix.py +++ b/pytential/symbolic/matrix.py @@ -41,6 +41,16 @@ def is_zero(x): def _resample_arg(queue, source, x): + """ + :arg queue: a :class:`pyopencl.CommandQueue`. + :arg source: a :class:`pytential.source.LayerPotentialSourceBase` subclass. + If it is not a layer potential source, no resampling is done. + :arg x: a :class:`numpy.ndarray`. + + :return: a resampled :class:`numpy.ndarray` (see + :method:`pytential.source.LayerPotentialSourceBase.resampler`). + """ + from pytential.source import LayerPotentialSourceBase if not isinstance(source, LayerPotentialSourceBase): return x @@ -59,6 +69,14 @@ def _resample_arg(queue, source, x): def _get_layer_potential_args(mapper, expr, source): + """ + :arg mapper: a :class:`pytential.symbolic.primitives.EvaluationMapperBase`. + :arg expr: symbolic layer potential expression. + :arg source: a :class:`pytential.source.LayerPotentialSourceBase`. + + :return: a mapping of kernel arguments *expr.kernel_arguments*. + """ + kernel_args = {} for arg_name, arg_expr in six.iteritems(expr.kernel_arguments): rec_arg = mapper.rec(arg_expr) @@ -68,6 +86,15 @@ def _get_layer_potential_args(mapper, expr, source): def _get_kernel_args(mapper, kernel, expr, source): + """ + :arg mapper: a :class:`pytential.symbolic.primitives.EvaluationMapperBase`. + :arg kernel: a :class:`sumpy.kernel.Kernel`. + :arg expr: symbolic layer potential expression. + :arg source: a :class:`pytential.source.LayerPotentialSourceBase`. + + :return: a mapping of kernel parameters included in *expr*. + """ + # NOTE: copied from pytential.symbolic.primitives.IntG inner_kernel_args = kernel.get_args() + kernel.get_source_args() inner_kernel_args = set(arg.loopy_arg.name for arg in inner_kernel_args) @@ -84,6 +111,14 @@ def _get_kernel_args(mapper, kernel, expr, source): def _get_weights_and_area_elements(queue, source, source_discr): + """ + :arg queue: a :class:`pyopencl.CommandQueue`. + :arg source: a :class:`pytential.source.LayerPotentialSourceBase` subclass. + :arg source_discr: a :class:`meshmode.discretization.Discretization`. + + :return: quadrature weights for each node in *source_discr* + """ + if source.quad_stage2_density_discr is source_discr: waa = source.weights_and_area_elements().with_queue(queue) else: @@ -98,6 +133,15 @@ def _get_weights_and_area_elements(queue, source, source_discr): def _get_centers_and_expansion_radii(queue, source, target_discr, qbx_forced_limit): + """ + :arg queue: a :class:`pyopencl.CommandQueue`. + :arg source: a :class:`pytential.source.LayerPotentialSourceBase` subclass. + :arg target_discr: a :class:`meshmode.discretization.Discretization`. + :arg qbx_forced_limit: an integer (*+1* or *-1*). + + :return: a tuple of `(centers, radii)` for each node in *target_discr*. + """ + if source.density_discr is target_discr: # NOTE: skip expensive target association from pytential.qbx.utils import get_centers_on_side -- GitLab From 96fb8e7e1d0f06521435ff0496779e9294792b3a Mon Sep 17 00:00:00 2001 From: Alex Fikl Date: Tue, 14 Aug 2018 20:53:28 -0500 Subject: [PATCH 24/34] matrix: fix issues with block matrix builders. Mainly added a more complicated operator in test_matrix. This exposed: * issues with the different map_variable implementations. * issues with resampling when unnneded, etc. --- pytential/symbolic/matrix.py | 276 ++++++++++++++++------------------- test/test_matrix.py | 43 +++--- 2 files changed, 153 insertions(+), 166 deletions(-) diff --git a/pytential/symbolic/matrix.py b/pytential/symbolic/matrix.py index 16396f6c..4a93421d 100644 --- a/pytential/symbolic/matrix.py +++ b/pytential/symbolic/matrix.py @@ -77,6 +77,13 @@ def _get_layer_potential_args(mapper, expr, source): :return: a mapping of kernel arguments *expr.kernel_arguments*. """ + # skip resampling if source and target are the same + from pytential.symbolic.primitives import DEFAULT_SOURCE, DEFAULT_TARGET + if ((expr.source is not DEFAULT_SOURCE) and + (expr.target is not DEFAULT_TARGET) and + (type(expr.source) is type(expr.target))): + source = None + kernel_args = {} for arg_name, arg_expr in six.iteritems(expr.kernel_arguments): rec_arg = mapper.rec(arg_expr) @@ -174,16 +181,12 @@ def _get_centers_and_expansion_radii(queue, source, target_discr, qbx_forced_lim # }}} +# {{{ base class for matrix builders -# {{{ QBX layer potential matrix builder - -# FIXME: PyOpenCL doesn't do all the required matrix math yet. -# We'll cheat and build the matrix on the host. - -class MatrixBuilder(EvaluationMapperBase): +class MatrixBuilderBase(EvaluationMapperBase): def __init__(self, queue, dep_expr, other_dep_exprs, dep_source, dep_discr, places, context): - super(MatrixBuilder, self).__init__(context=context) + super(MatrixBuilderBase, self).__init__(context=context) self.queue = queue self.dep_expr = dep_expr @@ -192,21 +195,38 @@ class MatrixBuilder(EvaluationMapperBase): self.dep_discr = dep_discr self.places = places + self.dep_nnodes = dep_discr.nnodes + + # {{{ + + def get_dep_variable(self): + return np.eye(self.dep_nnodes, dtype=np.float64) + + def is_kind_vector(self, x): + return len(x.shape) == 1 + + def is_kind_matrix(self, x): + return len(x.shape) == 2 + + # }}} + + # {{{ map_xxx implementation + def map_variable(self, expr): if expr == self.dep_expr: - return np.eye(self.dep_discr.nnodes, dtype=np.float64) + return self.get_dep_variable() elif expr in self.other_dep_exprs: return 0 else: - return super(MatrixBuilder, self).map_variable(expr) + return super(MatrixBuilderBase, self).map_variable(expr) def map_subscript(self, expr): if expr == self.dep_expr: - return np.eye(self.dep_discr.nnodes, dtype=np.float64) + return self.get_dep_variable() elif expr in self.other_dep_exprs: return 0 else: - return super(MatrixBuilder, self).map_subscript(expr) + return super(MatrixBuilderBase, self).map_subscript(expr) def map_sum(self, expr): sum_kind = None @@ -223,13 +243,12 @@ class MatrixBuilder(EvaluationMapperBase): continue if isinstance(rec_child, np.ndarray): - if len(rec_child.shape) == 2: + if self.is_kind_matrix(rec_child): term_kind = term_kind_matrix - elif len(rec_child.shape) == 1: + elif self.is_kind_vector(rec_child): term_kind = term_kind_vector else: raise RuntimeError("unexpected array rank") - else: term_kind = term_kind_scalar @@ -248,34 +267,84 @@ class MatrixBuilder(EvaluationMapperBase): mat_result = None vecs_and_scalars = 1 - for term in expr.children: - rec_term = self.rec(term) + for child in expr.children: + rec_child = self.rec(child) - if is_zero(rec_term): + if is_zero(rec_child): return 0 - if isinstance(rec_term, (np.number, int, float, complex)): - vecs_and_scalars = vecs_and_scalars * rec_term - elif isinstance(rec_term, np.ndarray): - if len(rec_term.shape) == 2: + if isinstance(rec_child, (np.number, int, float, complex)): + vecs_and_scalars = vecs_and_scalars * rec_child + elif isinstance(rec_child, np.ndarray): + if self.is_kind_matrix(rec_child): if mat_result is not None: raise RuntimeError("expression is nonlinear in %s" % self.dep_expr) else: - mat_result = rec_term + mat_result = rec_child else: - vecs_and_scalars = vecs_and_scalars * rec_term + vecs_and_scalars = vecs_and_scalars * rec_child if mat_result is not None: - if ( - isinstance(vecs_and_scalars, np.ndarray) - and len(vecs_and_scalars.shape) == 1): + if (isinstance(vecs_and_scalars, np.ndarray) + and self.is_kind_vector(vecs_and_scalars)): vecs_and_scalars = vecs_and_scalars[:, np.newaxis] return mat_result * vecs_and_scalars else: return vecs_and_scalars + def map_num_reference_derivative(self, expr): + rec_operand = self.rec(expr.operand) + + assert isinstance(rec_operand, np.ndarray) + if self.is_kind_matrix(rec_operand): + raise NotImplementedError("derivatives") + + where_discr = self.places[expr.where] + op = sym.NumReferenceDerivative(expr.ref_axes, sym.var("u")) + return bind(where_discr, op)( + self.queue, u=cl.array.to_device(self.queue, rec_operand)).get() + + def map_node_coordinate_component(self, expr): + where_discr = self.places[expr.where] + op = sym.NodeCoordinateComponent(expr.ambient_axis) + return bind(where_discr, op)(self.queue).get() + + def map_call(self, expr): + arg, = expr.parameters + rec_arg = self.rec(arg) + + if isinstance(rec_arg, np.ndarray) and self.is_kind_matrix(rec_arg): + raise RuntimeError("expression is nonlinear in variable") + + if isinstance(rec_arg, np.ndarray): + rec_arg = cl.array.to_device(self.queue, rec_arg) + + op = expr.function(sym.var("u")) + result = bind(self.dep_source, op)(self.queue, u=rec_arg) + + if isinstance(result, cl.array.Array): + result = result.get() + + return result + + # }}} + +# }}} + + +# {{{ QBX layer potential matrix builder + +# FIXME: PyOpenCL doesn't do all the required matrix math yet. +# We'll cheat and build the matrix on the host. + +class MatrixBuilder(MatrixBuilderBase): + def __init__(self, queue, dep_expr, other_dep_exprs, + dep_source, dep_discr, places, context): + super(MatrixBuilder, self).__init__(queue, dep_expr, other_dep_exprs, + dep_source, dep_discr, places, context) + def map_int_g(self, expr): # TODO: should this go into QBXPreprocessor / LocationTagger? where_source = expr.source @@ -291,7 +360,7 @@ class MatrixBuilder(EvaluationMapperBase): return 0 assert isinstance(rec_density, np.ndarray) - if len(rec_density.shape) != 2: + if not self.is_kind_matrix(rec_density): raise NotImplementedError("layer potentials on non-variables") kernel = expr.kernel @@ -320,6 +389,7 @@ class MatrixBuilder(EvaluationMapperBase): mat[:, :] *= waa.get(self.queue) if target_discr.nnodes != source_discr.nnodes: + # NOTE: we only resample sources assert target_discr.nnodes < source_discr.nnodes resampler = source.direct_resampler @@ -330,49 +400,12 @@ class MatrixBuilder(EvaluationMapperBase): return mat - # IntGdSource should have been removed by a preprocessor - - def map_num_reference_derivative(self, expr): - rec_operand = self.rec(expr.operand) - - assert isinstance(rec_operand, np.ndarray) - if len(rec_operand.shape) == 2: - raise NotImplementedError("derivatives") - - where_discr = self.places[expr.where] - op = sym.NumReferenceDerivative(expr.ref_axes, sym.var("u")) - return bind(where_discr, op)( - self.queue, u=cl.array.to_device(self.queue, rec_operand)).get() - - def map_node_coordinate_component(self, expr): - where_discr = self.places[expr.where] - op = sym.NodeCoordinateComponent(expr.ambient_axis) - return bind(where_discr, op)(self.queue).get() - - def map_call(self, expr): - arg, = expr.parameters - rec_arg = self.rec(arg) - - if isinstance(rec_arg, np.ndarray) and len(rec_arg.shape) == 2: - raise RuntimeError("expression is nonlinear in variable") - - if isinstance(rec_arg, np.ndarray): - rec_arg = cl.array.to_device(self.queue, rec_arg) - - op = expr.function(sym.var("u")) - result = bind(self.dep_source, op)(self.queue, u=rec_arg) - - if isinstance(result, cl.array.Array): - result = result.get() - - return result - # }}} # {{{ p2p matrix builder -class P2PMatrixBuilder(MatrixBuilder): +class P2PMatrixBuilder(MatrixBuilderBase): def __init__(self, queue, dep_expr, other_dep_exprs, dep_source, dep_discr, places, context, exclude_self=True): super(P2PMatrixBuilder, self).__init__(queue, @@ -391,7 +424,7 @@ class P2PMatrixBuilder(MatrixBuilder): return 0 assert isinstance(rec_density, np.ndarray) - if len(rec_density.shape) != 2: + if not self.is_kind_matrix(rec_density): raise NotImplementedError("layer potentials on non-variables") kernel = expr.kernel.get_base_kernel() @@ -418,72 +451,33 @@ class P2PMatrixBuilder(MatrixBuilder): # {{{ block matrix builders -class MatrixBlockBuilderBase(EvaluationMapperBase): +class MatrixBlockBuilderBase(MatrixBuilderBase): def __init__(self, queue, dep_expr, other_dep_exprs, dep_source, dep_discr, places, index_set, context): - super(MatrixBlockBuilderBase, self).__init__(context=context) + super(MatrixBlockBuilderBase, self).__init__(queue, + dep_expr, other_dep_exprs, dep_source, dep_discr, + places, context) - self.queue = queue - self.dep_expr = dep_expr - self.other_dep_exprs = other_dep_exprs - self.dep_source = dep_source - self.dep_discr = dep_discr - self.places = places self.index_set = index_set + self.dep_nnodes = index_set.col.indices.size - def _map_dep_variable(self): - return np.eye(self.index_set.col.indices.shape[0]) - - def map_variable(self, expr): - if expr == self.dep_expr: - return self._map_dep_variable() - elif expr in self.other_dep_exprs: - return 0 - else: - return super(MatrixBlockBuilderBase, self).map_variable(expr) - - def map_subscript(self, expr): - if expr == self.dep_expr: - return self.variable_identity() - elif expr in self.other_dep_exprs: - return 0 - else: - return super(MatrixBlockBuilderBase, self).map_subscript(expr) - - def map_num_reference_derivative(self, expr): - rec_operand = self.rec(expr.operand) - - assert isinstance(rec_operand, np.ndarray) - if len(rec_operand.shape) == 2: - raise NotImplementedError("derivatives") - - where_discr = self.places[expr.where] - op = sym.NumReferenceDerivative(expr.ref_axes, sym.var("u")) - return bind(where_discr, op)( - self.queue, u=cl.array.to_device(self.queue, rec_operand)).get() - - def map_node_coordinate_component(self, expr): - where_discr = self.places[expr.where] - op = sym.NodeCoordinateComponent(expr.ambient_axis) - return bind(where_discr, op)(self.queue).get() - - def map_call(self, expr): - arg, = expr.parameters - rec_arg = self.rec(arg) - - if isinstance(rec_arg, np.ndarray) and len(rec_arg.shape) == 2: - raise RuntimeError("expression is nonlinear in variable") - - if isinstance(rec_arg, np.ndarray): - rec_arg = cl.array.to_device(self.queue, rec_arg) + def get_dep_variable(self): + # NOTE: blocks are stored linearly, so the identity matrix for the + # variables themselves also needs to be flattened + tgtindices = self.index_set.linear_row_indices.get(self.queue) + srcindices = self.index_set.linear_col_indices.get(self.queue) - op = expr.function(sym.var("u")) - result = bind(self.dep_source, op)(self.queue, u=rec_arg) + return np.equal(tgtindices, srcindices).astype(np.float64) - if isinstance(result, cl.array.Array): - result = result.get() + def is_kind_vector(self, x): + # NOTE: since matrices are flattened, the only way to differentiate + # them from a vector is by size + return x.size == self.index_set.row.indices.size - return result + def is_kind_matrix(self, x): + # NOTE: since matrices are flattened, we recognize them by checking + # if they have the right size + return x.size == self.index_set.linear_row_indices.size class NearFieldBlockBuilder(MatrixBlockBuilderBase): @@ -493,15 +487,9 @@ class NearFieldBlockBuilder(MatrixBlockBuilderBase): dep_expr, other_dep_exprs, dep_source, dep_discr, places, index_set, context) - self.dummy = MatrixBlockBuilderBase(queue, + self.dummy = MatrixBuilderBase(queue, dep_expr, other_dep_exprs, dep_source, dep_discr, - places, index_set, context) - - def _map_dep_variable(self): - tgtindices = self.index_set.row.indices.get(self.queue).reshape(-1, 1) - srcindices = self.index_set.col.indices.get(self.queue).reshape(1, -1) - - return np.equal(tgtindices, srcindices).astype(np.float64) + places, context) def map_int_g(self, expr): source = self.places[expr.source] @@ -511,16 +499,16 @@ class NearFieldBlockBuilder(MatrixBlockBuilderBase): if source_discr is not target_discr: raise NotImplementedError() - rec_density = self.dummy.rec(expr.density) + rec_density = self.rec(expr.density) if is_zero(rec_density): return 0 assert isinstance(rec_density, np.ndarray) - if len(rec_density.shape) != 2: + if not self.is_kind_matrix(rec_density): raise NotImplementedError("layer potentials on non-variables") kernel = expr.kernel - kernel_args = _get_layer_potential_args(self, expr, source) + kernel_args = _get_layer_potential_args(self.dummy, expr, source) from sumpy.expansion.local import LineTaylorLocalExpansion local_expn = LineTaylorLocalExpansion(kernel, source.qbx_order) @@ -557,32 +545,26 @@ class FarFieldBlockBuilder(MatrixBlockBuilderBase): dep_expr, other_dep_exprs, dep_source, dep_discr, places, index_set, context) - self.dummy = MatrixBlockBuilderBase(queue, - dep_expr, other_dep_exprs, dep_source, dep_discr, - places, index_set, context) self.exclude_self = exclude_self - - def _map_dep_variable(self): - tgtindices = self.index_set.row.indices.get(self.queue).reshape(-1, 1) - srcindices = self.index_set.col.indices.get(self.queue).reshape(1, -1) - - return np.equal(tgtindices, srcindices).astype(np.float64) + self.dummy = MatrixBuilderBase(queue, + dep_expr, other_dep_exprs, dep_source, dep_discr, + places, context) def map_int_g(self, expr): source = self.places[expr.source] source_discr = self.places.get_discretization(expr.source) target_discr = self.places.get_discretization(expr.target) - rec_density = self.dummy.rec(expr.density) + rec_density = self.rec(expr.density) if is_zero(rec_density): return 0 assert isinstance(rec_density, np.ndarray) - if len(rec_density.shape) != 2: + if not self.is_kind_matrix(rec_density): raise NotImplementedError("layer potentials on non-variables") kernel = expr.kernel.get_base_kernel() - kernel_args = _get_kernel_args(self, kernel, expr, source) + kernel_args = _get_kernel_args(self.dummy, kernel, expr, source) if self.exclude_self: kernel_args["target_to_source"] = \ cl.array.arange(self.queue, 0, target_discr.nnodes, dtype=np.int) diff --git a/test/test_matrix.py b/test/test_matrix.py index d16dd32a..5aaae91f 100644 --- a/test/test_matrix.py +++ b/test/test_matrix.py @@ -50,7 +50,7 @@ from pyopencl.tools import ( # noqa def _build_op(lpot_id, k=0, ndim=2, - qbx_forced_limit=-1): + qbx_forced_limit="avg"): from sumpy.kernel import LaplaceKernel, HelmholtzKernel if k: @@ -76,15 +76,16 @@ def _build_op(lpot_id, u0_sym, u1_sym = u_sym op = make_obj_array([ - sym.Sp(knl, u0_sym, **lpot_kwargs) + - sym.D(knl, u1_sym, **lpot_kwargs), - - sym.S(knl, 0.4 * u0_sym, **lpot_kwargs) + - 0.3 * sym.D(knl, u0_sym, **lpot_kwargs) + sym.Sp(knl, u0_sym, **lpot_kwargs) + + sym.D(knl, u1_sym, **lpot_kwargs), + sym.S(knl, 0.4 * u0_sym, **lpot_kwargs) + + 0.3 * sym.D(knl, u0_sym, **lpot_kwargs) ]) else: raise ValueError("Unknown lpot_id: {}".format(lpot_id)) + op = 0.5 * u_sym + op + return op, u_sym, knl_kwargs @@ -94,7 +95,7 @@ def _build_op(lpot_id, @pytest.mark.parametrize("curve_f", [ partial(ellipse, 3), NArmedStarfish(5, 0.25)]) -@pytest.mark.parametrize("lpot_id", [1, 3]) +@pytest.mark.parametrize("lpot_id", [2, 3]) def test_matrix_build(ctx_factory, k, curve_f, lpot_id, visualize=False): cl_ctx = ctx_factory() queue = cl.CommandQueue(cl_ctx) @@ -212,26 +213,26 @@ def test_p2p_block_builder(ctx_factory, factor, ndim, lpot_id, from sumpy.tools import MatrixBlockIndexRanges index_set = MatrixBlockIndexRanges(ctx, tgtindices, srcindices) - from pytential.symbolic.matrix import FarFieldBlockBuilder - mbuilder = FarFieldBlockBuilder(queue, + from pytential.symbolic.matrix import P2PMatrixBuilder + mbuilder = P2PMatrixBuilder(queue, dep_expr=u_sym, other_dep_exprs=[], dep_source=places[domains[0]], dep_discr=places.get_discretization(domains[0]), places=places, - index_set=index_set, context={}) - blk = mbuilder(expr) + mat = mbuilder(expr) - from pytential.symbolic.matrix import P2PMatrixBuilder - mbuilder = P2PMatrixBuilder(queue, + from pytential.symbolic.matrix import FarFieldBlockBuilder + mbuilder = FarFieldBlockBuilder(queue, dep_expr=u_sym, other_dep_exprs=[], dep_source=places[domains[0]], dep_discr=places.get_discretization(domains[0]), places=places, + index_set=index_set, context={}) - mat = mbuilder(expr) + blk = mbuilder(expr) index_set = index_set.get(queue) if visualize and ndim == 2: @@ -263,9 +264,11 @@ def test_p2p_block_builder(ctx_factory, factor, ndim, lpot_id, assert error < eps +@pytest.mark.parametrize("factor", [1.0, 0.6]) @pytest.mark.parametrize("ndim", [2, 3]) -@pytest.mark.parametrize("lpot_id", [1]) -def test_qbx_block_builder(ctx_factory, ndim, lpot_id, visualize=False): +@pytest.mark.parametrize("lpot_id", [1, 2]) +def test_qbx_block_builder(ctx_factory, factor, ndim, lpot_id, + visualize=False): ctx = ctx_factory() queue = cl.CommandQueue(ctx) @@ -278,9 +281,6 @@ def test_qbx_block_builder(ctx_factory, ndim, lpot_id, visualize=False): qbx = _build_qbx_discr(queue, target_order=target_order, ndim=ndim) op, u_sym, _ = _build_op(lpot_id, ndim=ndim) - tgtindices = _build_block_index(qbx.density_discr) - srcindices = _build_block_index(qbx.density_discr) - # NOTE: NearFieldBlockBuilder only does stage1/stage1 or stage2/stage2, # so we need to hardcode the discr for MatrixBuilder too, since the # defaults are different @@ -289,8 +289,13 @@ def test_qbx_block_builder(ctx_factory, ndim, lpot_id, visualize=False): from pytential.symbolic.execution import GeometryCollection, _prepare_expr places = GeometryCollection(qbx, auto_where=where) expr = _prepare_expr(places, op) + density_discr = places.get_discretization(where[0]) from sumpy.tools import MatrixBlockIndexRanges + srcindices = _build_block_index(density_discr, + method='nodes', factor=factor) + tgtindices = _build_block_index(density_discr, + method='nodes', factor=factor) index_set = MatrixBlockIndexRanges(ctx, tgtindices, srcindices) from pytential.symbolic.matrix import NearFieldBlockBuilder -- GitLab From 836d3d0f6d46219ceeb322f5e805825376d39caa Mon Sep 17 00:00:00 2001 From: Alex Fikl Date: Tue, 14 Aug 2018 21:04:49 -0500 Subject: [PATCH 25/34] matrix: flake8 fixes --- pytential/symbolic/execution.py | 20 ++++++++++---------- pytential/symbolic/matrix.py | 7 ++++--- test/test_matrix.py | 8 ++++---- 3 files changed, 18 insertions(+), 17 deletions(-) diff --git a/pytential/symbolic/execution.py b/pytential/symbolic/execution.py index b1e8eb4a..4ebd2831 100644 --- a/pytential/symbolic/execution.py +++ b/pytential/symbolic/execution.py @@ -498,16 +498,16 @@ def _bmat(blocks, dtypes): ncolumns = blocks.shape[1] # "block row starts"/"block column starts" - brs = np.cumsum([0] + - [single_valued(blocks[ibrow, ibcol].shape[0] - for ibcol in range(ncolumns) - if not is_zero(blocks[ibrow, ibcol])) + brs = np.cumsum([0] + + [single_valued(blocks[ibrow, ibcol].shape[0] + for ibcol in range(ncolumns) + if not is_zero(blocks[ibrow, ibcol])) for ibrow in range(nrows)]) - bcs = np.cumsum([0] + - [single_valued(blocks[ibrow, ibcol].shape[1] - for ibrow in range(nrows) - if not is_zero(blocks[ibrow, ibcol])) + bcs = np.cumsum([0] + + [single_valued(blocks[ibrow, ibcol].shape[1] + for ibrow in range(nrows) + if not is_zero(blocks[ibrow, ibcol])) for ibcol in range(ncolumns)]) result = np.zeros((brs[-1], bcs[-1]), @@ -571,8 +571,8 @@ def build_matrix(queue, places, exprs, input_exprs, domains=None, mbuilder = MatrixBuilder( queue, dep_expr=input_exprs[ibcol], - other_dep_exprs=(input_exprs[:ibcol] + - input_exprs[ibcol + 1:]), + other_dep_exprs=(input_exprs[:ibcol] + + input_exprs[ibcol + 1:]), dep_source=places[domains[ibcol]], dep_discr=places.get_discretization(domains[ibcol]), places=places, diff --git a/pytential/symbolic/matrix.py b/pytential/symbolic/matrix.py index 4a93421d..6a1e619b 100644 --- a/pytential/symbolic/matrix.py +++ b/pytential/symbolic/matrix.py @@ -79,9 +79,9 @@ def _get_layer_potential_args(mapper, expr, source): # skip resampling if source and target are the same from pytential.symbolic.primitives import DEFAULT_SOURCE, DEFAULT_TARGET - if ((expr.source is not DEFAULT_SOURCE) and - (expr.target is not DEFAULT_TARGET) and - (type(expr.source) is type(expr.target))): + if ((expr.source is not DEFAULT_SOURCE) + and (expr.target is not DEFAULT_TARGET) + and (isinstance(expr.source, type(expr.target)))): source = None kernel_args = {} @@ -181,6 +181,7 @@ def _get_centers_and_expansion_radii(queue, source, target_discr, qbx_forced_lim # }}} + # {{{ base class for matrix builders class MatrixBuilderBase(EvaluationMapperBase): diff --git a/test/test_matrix.py b/test/test_matrix.py index 5aaae91f..0c26019c 100644 --- a/test/test_matrix.py +++ b/test/test_matrix.py @@ -256,8 +256,8 @@ def test_p2p_block_builder(ctx_factory, factor, ndim, lpot_id, for i in range(index_set.nblocks): eps = 1.0e-14 * la.norm(index_set.take(mat, i)) - error = la.norm(index_set.block_take(blk, i) - - index_set.take(mat, i)) + error = la.norm(index_set.block_take(blk, i) + - index_set.take(mat, i)) if visualize: print('block[{:04}]: {:.5e}'.format(i, error)) @@ -341,8 +341,8 @@ def test_qbx_block_builder(ctx_factory, factor, ndim, lpot_id, for i in range(index_set.nblocks): eps = 1.0e-14 * la.norm(index_set.take(mat, i)) - error = la.norm(index_set.block_take(blk, i) - - index_set.take(mat, i)) + error = la.norm(index_set.block_take(blk, i) + - index_set.take(mat, i)) if visualize: print('block[{:04}]: {:.5e}'.format(i, error)) -- GitLab From 7f1416a3e72172d2611b5b2073892528a00b2caf Mon Sep 17 00:00:00 2001 From: Alex Fikl Date: Wed, 15 Aug 2018 20:05:29 -0500 Subject: [PATCH 26/34] execution: add some docs to GeometryCollection --- pytential/symbolic/execution.py | 76 ++++++++++++++++++++++++--------- pytential/symbolic/matrix.py | 59 +++++++++++++------------ test/test_matrix.py | 12 +++--- 3 files changed, 96 insertions(+), 51 deletions(-) diff --git a/pytential/symbolic/execution.py b/pytential/symbolic/execution.py index 4ebd2831..6ca61346 100644 --- a/pytential/symbolic/execution.py +++ b/pytential/symbolic/execution.py @@ -286,7 +286,7 @@ class MatVecOp: def _prepare_domains(nresults, places, domains, default_domain): """ :arg nresults: number of results. - :arg places: :class:`pytential.symbolic.execution.GeometryCollection`. + :arg places: a :class:`pytential.symbolic.execution.GeometryCollection`. :arg domains: recommended domains. :arg default_domain: default value for domains which are not provided. @@ -302,7 +302,7 @@ def _prepare_domains(nresults, places, domains, default_domain): "default domain to be defined in places") dom_name = default_domain return nresults * [dom_name] - if not isinstance(domains, (list, tuple)): + elif not isinstance(domains, (list, tuple)): dom_name = domains return nresults * [dom_name] @@ -337,10 +337,34 @@ def _prepare_expr(places, expr): # {{{ bound expression class GeometryCollection(object): + """A collection of geometry-related objects. This class is meant to hold + a specific combination of sources and targets and cache any common + expressions specific to it, e.g. metric terms, etc. + + .. method:: __getitem__ + .. method:: get_discretization + .. method:: get_cache + """ + def __init__(self, places, auto_where=None): + """ + :arg places: a scalar, tuple of or mapping of symbolic names to + geometry objects. Supported objects are + :class:`~pytential.source.PotentialSource`, + :class:`~potential.target.TargetBase` and + :class:`~meshmode.discretization.Discretization`. + :arg auto_where: location identifier for each geometry object, used + to denote specific discretizations, e.g. in the case where + *places* is a :class:`~pytential.source.LayerPotentialSourceBase`. + By default, we assume + :class:`~pytential.symbolic.primitives.DEFAULT_SOURCE` and + :class:`~pytential.symbolic.primitives.DEFAULT_TARGET` for + sources and targets, respectively. + """ + + from pytential.target import TargetBase from meshmode.discretization import Discretization from pytential.source import LayerPotentialSourceBase, PotentialSource - from pytential.target import TargetBase if auto_where is None: auto_where = DEFAULT_SOURCE, DEFAULT_TARGET @@ -350,7 +374,7 @@ class GeometryCollection(object): if isinstance(places, LayerPotentialSourceBase): self.places[self.source] = places self.places[self.target] = \ - self.get_discretization(self.target, lpot=places) + self._get_lpot_discretization(self.target, places) elif isinstance(places, (Discretization, TargetBase)): self.places[self.target] = places elif isinstance(places, tuple): @@ -375,20 +399,14 @@ class GeometryCollection(object): def target(self): return self.where[1] - def get_discretization(self, where, lpot=None): - if where is DEFAULT_SOURCE or where is DEFAULT_TARGET: - where = QBXSourceStage1(where) - - if lpot is None: - if where in self.places: - lpot = self.places[where] - else: - lpot = self.places.get(where.where, None) - + def _get_lpot_discretization(self, where, lpot): from pytential.source import LayerPotentialSourceBase if not isinstance(lpot, LayerPotentialSourceBase): return lpot + if where is DEFAULT_SOURCE or where is DEFAULT_TARGET: + where = QBXSourceStage1(where) + if isinstance(where, QBXSourceStage1): return lpot.density_discr if isinstance(where, QBXSourceStage2): @@ -396,7 +414,28 @@ class GeometryCollection(object): if isinstance(where, QBXSourceQuadStage2): return lpot.quad_stage2_density_discr - raise ValueError('Unknown `where` identifier: {}'.format(type(where))) + raise ValueError('unknown `where` identifier: {}'.format(type(where))) + + def get_discretization(self, where): + """ + :arg where: location identifier. + + :return: a geometry object in the collection corresponding to the + key *where*. If it is a + :class:`~pytential.source.LayerPotentialSourceBase`, we look for + the corresponding :class:`~meshmode.discretization.Discretization` + in its attributes instead. + """ + + if where in self.places: + lpot = self.places[where] + else: + lpot = self.places.get(getattr(where, 'where', None), None) + + if lpot is None: + raise KeyError('`where` not in the collection: {}'.format(where)) + + return self._get_lpot_discretization(where, lpot) def __getitem__(self, where): return self.places[where] @@ -471,10 +510,9 @@ class BoundExpression(object): def bind(places, expr, auto_where=None): """ - :arg places: a collection or mapping of symbolic names to - :class:`meshmode.discretization.Discretization` objects, subclasses - of :class:`pytential.source.PotentialSource` or subclasses of - :class:`pytential.target.TargetBase`. + :arg places: a :class:`pytential.symbolic.execution.GeometryCollection`. + Alternatively, any list or mapping that is a valid argument for its + constructor can also be used. :arg auto_where: for simple source-to-self or source-to-target evaluations, find 'where' attributes automatically. """ diff --git a/pytential/symbolic/matrix.py b/pytential/symbolic/matrix.py index 6a1e619b..01e60caf 100644 --- a/pytential/symbolic/matrix.py +++ b/pytential/symbolic/matrix.py @@ -70,11 +70,11 @@ def _resample_arg(queue, source, x): def _get_layer_potential_args(mapper, expr, source): """ - :arg mapper: a :class:`pytential.symbolic.primitives.EvaluationMapperBase`. + :arg mapper: a :class:`pytential.symbolic.matrix.MatrixBuilderBase`. :arg expr: symbolic layer potential expression. :arg source: a :class:`pytential.source.LayerPotentialSourceBase`. - :return: a mapping of kernel arguments *expr.kernel_arguments*. + :return: a mapping of kernel arguments evaluated by the *mapper*. """ # skip resampling if source and target are the same @@ -94,12 +94,12 @@ def _get_layer_potential_args(mapper, expr, source): def _get_kernel_args(mapper, kernel, expr, source): """ - :arg mapper: a :class:`pytential.symbolic.primitives.EvaluationMapperBase`. + :arg mapper: a :class:`pytential.symbolic.matrix.MatrixBuilderBase`. :arg kernel: a :class:`sumpy.kernel.Kernel`. :arg expr: symbolic layer potential expression. :arg source: a :class:`pytential.source.LayerPotentialSourceBase`. - :return: a mapping of kernel parameters included in *expr*. + :return: a mapping of kernel arguments evaluated by the *mapper*. """ # NOTE: copied from pytential.symbolic.primitives.IntG @@ -120,10 +120,10 @@ def _get_kernel_args(mapper, kernel, expr, source): def _get_weights_and_area_elements(queue, source, source_discr): """ :arg queue: a :class:`pyopencl.CommandQueue`. - :arg source: a :class:`pytential.source.LayerPotentialSourceBase` subclass. + :arg source: a :class:`pytential.source.LayerPotentialSourceBase`. :arg source_discr: a :class:`meshmode.discretization.Discretization`. - :return: quadrature weights for each node in *source_discr* + :return: quadrature weights for each node in *source_discr*. """ if source.quad_stage2_density_discr is source_discr: @@ -142,7 +142,7 @@ def _get_weights_and_area_elements(queue, source, source_discr): def _get_centers_and_expansion_radii(queue, source, target_discr, qbx_forced_limit): """ :arg queue: a :class:`pyopencl.CommandQueue`. - :arg source: a :class:`pytential.source.LayerPotentialSourceBase` subclass. + :arg source: a :class:`pytential.source.LayerPotentialSourceBase`. :arg target_discr: a :class:`meshmode.discretization.Discretization`. :arg qbx_forced_limit: an integer (*+1* or *-1*). @@ -187,6 +187,20 @@ def _get_centers_and_expansion_radii(queue, source, target_discr, qbx_forced_lim class MatrixBuilderBase(EvaluationMapperBase): def __init__(self, queue, dep_expr, other_dep_exprs, dep_source, dep_discr, places, context): + """ + :arg queue: a :class:`pyopencl.CommandQueue`. + :arg dep_expr: symbolic expression for the input block column + that the builder is evaluating. + :arg other_dep_exprs: symbolic expressions for the remaining input + block columns. + :arg dep_source: a :class:`pytential.source.LayerPotentialSourceBase` + for the given *dep_expr*. + :arg dep_discr: a concerete :class:`meshmode.discretization.Discretization` + for the given *dep_expr*. + :arg places: a :class:`pytential.symbolic.execution.GeometryCollection` + for all the sources and targets the builder is expected to + encounter. + """ super(MatrixBuilderBase, self).__init__(context=context) self.queue = queue @@ -347,7 +361,6 @@ class MatrixBuilder(MatrixBuilderBase): dep_source, dep_discr, places, context) def map_int_g(self, expr): - # TODO: should this go into QBXPreprocessor / LocationTagger? where_source = expr.source if where_source is sym.DEFAULT_SOURCE: where_source = sym.QBXSourceQuadStage2(expr.source) @@ -463,12 +476,10 @@ class MatrixBlockBuilderBase(MatrixBuilderBase): self.dep_nnodes = index_set.col.indices.size def get_dep_variable(self): - # NOTE: blocks are stored linearly, so the identity matrix for the - # variables themselves also needs to be flattened - tgtindices = self.index_set.linear_row_indices.get(self.queue) - srcindices = self.index_set.linear_col_indices.get(self.queue) - - return np.equal(tgtindices, srcindices).astype(np.float64) + # NOTE: block builders only support identity operators acting on + # the density. If other operators are needed, the user is meant to + # do the required matrix-matrix products + return 1.0 def is_kind_vector(self, x): # NOTE: since matrices are flattened, the only way to differentiate @@ -504,9 +515,8 @@ class NearFieldBlockBuilder(MatrixBlockBuilderBase): if is_zero(rec_density): return 0 - assert isinstance(rec_density, np.ndarray) - if not self.is_kind_matrix(rec_density): - raise NotImplementedError("layer potentials on non-variables") + if not np.isscalar(rec_density): + raise NotImplementedError() kernel = expr.kernel kernel_args = _get_layer_potential_args(self.dummy, expr, source) @@ -532,16 +542,14 @@ class NearFieldBlockBuilder(MatrixBlockBuilderBase): waa = _get_weights_and_area_elements(self.queue, source, source_discr) mat *= waa[self.index_set.linear_col_indices] - mat = mat.get(self.queue) - - # TODO: multiply with rec_density + mat = rec_density * mat.get(self.queue) return mat class FarFieldBlockBuilder(MatrixBlockBuilderBase): def __init__(self, queue, dep_expr, other_dep_exprs, dep_source, dep_discr, - places, index_set, context, exclude_self=True): + places, index_set, context, exclude_self=False): super(FarFieldBlockBuilder, self).__init__(queue, dep_expr, other_dep_exprs, dep_source, dep_discr, places, index_set, context) @@ -560,9 +568,8 @@ class FarFieldBlockBuilder(MatrixBlockBuilderBase): if is_zero(rec_density): return 0 - assert isinstance(rec_density, np.ndarray) - if not self.is_kind_matrix(rec_density): - raise NotImplementedError("layer potentials on non-variables") + if not np.isscalar(rec_density): + raise NotImplementedError() kernel = expr.kernel.get_base_kernel() kernel_args = _get_kernel_args(self.dummy, kernel, expr, source) @@ -579,9 +586,7 @@ class FarFieldBlockBuilder(MatrixBlockBuilderBase): sources=source_discr.nodes(), index_set=self.index_set, **kernel_args) - mat = mat.get() - - # TODO: need to multiply by rec_density + mat = rec_density * mat.get(self.queue) return mat diff --git a/test/test_matrix.py b/test/test_matrix.py index 0c26019c..69f923db 100644 --- a/test/test_matrix.py +++ b/test/test_matrix.py @@ -65,11 +65,11 @@ def _build_op(lpot_id, if lpot_id == 1: # scalar single-layer potential u_sym = sym.var("u") - op = sym.S(knl, u_sym, **lpot_kwargs) + op = sym.S(knl, 0.3 * u_sym, **lpot_kwargs) elif lpot_id == 2: # scalar double-layer potential u_sym = sym.var("u") - op = sym.D(knl, u_sym, **lpot_kwargs) + op = sym.D(knl, 0.3 * u_sym, **lpot_kwargs) elif lpot_id == 3: # vector potential u_sym = sym.make_sym_vector("u", 2) @@ -84,7 +84,7 @@ def _build_op(lpot_id, else: raise ValueError("Unknown lpot_id: {}".format(lpot_id)) - op = 0.5 * u_sym + op + # op = 0.5 * u_sym + op return op, u_sym, knl_kwargs @@ -220,7 +220,8 @@ def test_p2p_block_builder(ctx_factory, factor, ndim, lpot_id, dep_source=places[domains[0]], dep_discr=places.get_discretization(domains[0]), places=places, - context={}) + context={}, + exclude_self=True) mat = mbuilder(expr) from pytential.symbolic.matrix import FarFieldBlockBuilder @@ -231,7 +232,8 @@ def test_p2p_block_builder(ctx_factory, factor, ndim, lpot_id, dep_discr=places.get_discretization(domains[0]), places=places, index_set=index_set, - context={}) + context={}, + exclude_self=True) blk = mbuilder(expr) index_set = index_set.get(queue) -- GitLab From 4c6373744c8a23b59fe261ee9c3e244a3bd3a0ed Mon Sep 17 00:00:00 2001 From: Alex Fikl Date: Wed, 15 Aug 2018 20:31:50 -0500 Subject: [PATCH 27/34] matrix: make sure block builders don't allow composition --- pytential/symbolic/matrix.py | 121 +++++++++++++++++++++++------------ test/test_matrix.py | 9 +-- 2 files changed, 85 insertions(+), 45 deletions(-) diff --git a/pytential/symbolic/matrix.py b/pytential/symbolic/matrix.py index 01e60caf..77cc31ae 100644 --- a/pytential/symbolic/matrix.py +++ b/pytential/symbolic/matrix.py @@ -346,6 +346,44 @@ class MatrixBuilderBase(EvaluationMapperBase): # }}} + +class MatrixBlockBuilderBase(MatrixBuilderBase): + """Evaluate individual blocks of a matrix operator. + + Unlike, e.g. :class:`MatrixBuilder`, matrix block builders are + significantly reduced in scope. They are basically just meant + to evaluate linear combinations of layer potential operators. + For example, they do not support composition of operators because we + assume that each operator acts directly on the density. + """ + + def __init__(self, queue, dep_expr, other_dep_exprs, + dep_source, dep_discr, places, index_set, context): + """ + :arg index_set: a :class:`sumpy.tools.MatrixBlockIndexRanges` class + describing which blocks are going to be evaluated. + """ + + super(MatrixBlockBuilderBase, self).__init__(queue, + dep_expr, other_dep_exprs, dep_source, dep_discr, + places, context) + + self.index_set = index_set + self.dep_nnodes = index_set.col.indices.size + + def get_dep_variable(self): + return 1.0 + + def is_kind_vector(self, x): + # NOTE: since matrices are flattened, the only way to differentiate + # them from a vector is by size + return x.size == self.index_set.row.indices.size + + def is_kind_matrix(self, x): + # NOTE: since matrices are flattened, we recognize them by checking + # if they have the right size + return x.size == self.index_set.linear_row_indices.size + # }}} @@ -465,43 +503,31 @@ class P2PMatrixBuilder(MatrixBuilderBase): # {{{ block matrix builders -class MatrixBlockBuilderBase(MatrixBuilderBase): - def __init__(self, queue, dep_expr, other_dep_exprs, - dep_source, dep_discr, places, index_set, context): - super(MatrixBlockBuilderBase, self).__init__(queue, - dep_expr, other_dep_exprs, dep_source, dep_discr, - places, context) - - self.index_set = index_set - self.dep_nnodes = index_set.col.indices.size - - def get_dep_variable(self): - # NOTE: block builders only support identity operators acting on - # the density. If other operators are needed, the user is meant to - # do the required matrix-matrix products - return 1.0 - - def is_kind_vector(self, x): - # NOTE: since matrices are flattened, the only way to differentiate - # them from a vector is by size - return x.size == self.index_set.row.indices.size - - def is_kind_matrix(self, x): - # NOTE: since matrices are flattened, we recognize them by checking - # if they have the right size - return x.size == self.index_set.linear_row_indices.size - - class NearFieldBlockBuilder(MatrixBlockBuilderBase): def __init__(self, queue, dep_expr, other_dep_exprs, dep_source, dep_discr, places, index_set, context): super(NearFieldBlockBuilder, self).__init__(queue, - dep_expr, other_dep_exprs, dep_source, dep_discr, - places, index_set, context) + dep_expr, other_dep_exprs, dep_source, dep_discr, + places, index_set, context) + + # NOTE: we need additional mappers to redirect some operations: + # * mat_mapper is used to compute any kernel arguments that need to + # be computed on the full discretization, ignoring our index_set, + # e.g the normal in a double layer potential + # * blk_mapper is used to recursively compute the density to + # a layer potential operator to ensure there is no composition + self.mat_mapper = MatrixBuilderBase(queue, + dep_expr, other_dep_exprs, dep_source, dep_discr, + places, context) + self.blk_mapper = MatrixBlockBuilderBase(queue, + dep_expr, other_dep_exprs, dep_source, dep_discr, + places, index_set, context) - self.dummy = MatrixBuilderBase(queue, - dep_expr, other_dep_exprs, dep_source, dep_discr, - places, context) + def get_dep_variable(self): + tgtindices = self.index_set.linear_row_indices.get(self.queue) + srcindices = self.index_set.linear_col_indices.get(self.queue) + + return np.equal(tgtindices, srcindices).astype(np.float64) def map_int_g(self, expr): source = self.places[expr.source] @@ -511,7 +537,7 @@ class NearFieldBlockBuilder(MatrixBlockBuilderBase): if source_discr is not target_discr: raise NotImplementedError() - rec_density = self.rec(expr.density) + rec_density = self.blk_mapper.rec(expr.density) if is_zero(rec_density): return 0 @@ -519,7 +545,7 @@ class NearFieldBlockBuilder(MatrixBlockBuilderBase): raise NotImplementedError() kernel = expr.kernel - kernel_args = _get_layer_potential_args(self.dummy, expr, source) + kernel_args = _get_layer_potential_args(self.mat_mapper, expr, source) from sumpy.expansion.local import LineTaylorLocalExpansion local_expn = LineTaylorLocalExpansion(kernel, source.qbx_order) @@ -551,20 +577,33 @@ class FarFieldBlockBuilder(MatrixBlockBuilderBase): def __init__(self, queue, dep_expr, other_dep_exprs, dep_source, dep_discr, places, index_set, context, exclude_self=False): super(FarFieldBlockBuilder, self).__init__(queue, - dep_expr, other_dep_exprs, dep_source, dep_discr, - places, index_set, context) + dep_expr, other_dep_exprs, dep_source, dep_discr, + places, index_set, context) + # NOTE: same mapper issues as in the NearFieldBlockBuilder self.exclude_self = exclude_self - self.dummy = MatrixBuilderBase(queue, - dep_expr, other_dep_exprs, dep_source, dep_discr, - places, context) + self.mat_mapper = MatrixBuilderBase(queue, + dep_expr, other_dep_exprs, dep_source, dep_discr, + places, context) + self.blk_mapper = MatrixBlockBuilderBase(queue, + dep_expr, other_dep_exprs, dep_source, dep_discr, + places, index_set, context) + + def get_dep_variable(self): + tgtindices = self.index_set.linear_row_indices.get(self.queue) + srcindices = self.index_set.linear_col_indices.get(self.queue) + + return np.equal(tgtindices, srcindices).astype(np.float64) def map_int_g(self, expr): source = self.places[expr.source] source_discr = self.places.get_discretization(expr.source) target_discr = self.places.get_discretization(expr.target) - rec_density = self.rec(expr.density) + if source_discr is not target_discr: + raise NotImplementedError() + + rec_density = self.blk_mapper.rec(expr.density) if is_zero(rec_density): return 0 @@ -572,7 +611,7 @@ class FarFieldBlockBuilder(MatrixBlockBuilderBase): raise NotImplementedError() kernel = expr.kernel.get_base_kernel() - kernel_args = _get_kernel_args(self.dummy, kernel, expr, source) + kernel_args = _get_kernel_args(self.mat_mapper, kernel, expr, source) if self.exclude_self: kernel_args["target_to_source"] = \ cl.array.arange(self.queue, 0, target_discr.nnodes, dtype=np.int) diff --git a/test/test_matrix.py b/test/test_matrix.py index 69f923db..02eaca0d 100644 --- a/test/test_matrix.py +++ b/test/test_matrix.py @@ -65,11 +65,12 @@ def _build_op(lpot_id, if lpot_id == 1: # scalar single-layer potential u_sym = sym.var("u") - op = sym.S(knl, 0.3 * u_sym, **lpot_kwargs) + op = sym.S(knl, u_sym, **lpot_kwargs) elif lpot_id == 2: - # scalar double-layer potential + # scalar combination of layer potentials u_sym = sym.var("u") - op = sym.D(knl, 0.3 * u_sym, **lpot_kwargs) + op = sym.S(knl, 0.3 * u_sym, **lpot_kwargs) \ + + sym.D(knl, 0.5 * u_sym, **lpot_kwargs) elif lpot_id == 3: # vector potential u_sym = sym.make_sym_vector("u", 2) @@ -84,7 +85,7 @@ def _build_op(lpot_id, else: raise ValueError("Unknown lpot_id: {}".format(lpot_id)) - # op = 0.5 * u_sym + op + op = 0.5 * u_sym + op return op, u_sym, knl_kwargs -- GitLab From 79369904e7c71e09485404715625fc62d594becd Mon Sep 17 00:00:00 2001 From: Alex Fikl Date: Wed, 15 Aug 2018 21:08:56 -0500 Subject: [PATCH 28/34] execution: remove some properties from GeometryCollection --- pytential/symbolic/execution.py | 37 +++++++++++++-------------------- test/test_matrix.py | 5 +++-- 2 files changed, 18 insertions(+), 24 deletions(-) diff --git a/pytential/symbolic/execution.py b/pytential/symbolic/execution.py index 6ca61346..71300d83 100644 --- a/pytential/symbolic/execution.py +++ b/pytential/symbolic/execution.py @@ -367,20 +367,22 @@ class GeometryCollection(object): from pytential.source import LayerPotentialSourceBase, PotentialSource if auto_where is None: - auto_where = DEFAULT_SOURCE, DEFAULT_TARGET - self.where = auto_where + source_where, target_where = DEFAULT_SOURCE, DEFAULT_TARGET + else: + source_where, target_where = auto_where + self.where = (source_where, target_where) self.places = {} if isinstance(places, LayerPotentialSourceBase): - self.places[self.source] = places - self.places[self.target] = \ - self._get_lpot_discretization(self.target, places) + self.places[source_where] = places + self.places[target_where] = \ + self._get_lpot_discretization(target_where, places) elif isinstance(places, (Discretization, TargetBase)): - self.places[self.target] = places + self.places[target_where] = places elif isinstance(places, tuple): source_discr, target_discr = places - self.places[self.source] = source_discr - self.places[self.target] = target_discr + self.places[source_where] = source_discr + self.places[target_where] = target_discr else: self.places = places.copy() @@ -391,14 +393,6 @@ class GeometryCollection(object): self.caches = {} - @property - def source(self): - return self.where[0] - - @property - def target(self): - return self.where[1] - def _get_lpot_discretization(self, where, lpot): from pytential.source import LayerPotentialSourceBase if not isinstance(lpot, LayerPotentialSourceBase): @@ -562,11 +556,9 @@ def build_matrix(queue, places, exprs, input_exprs, domains=None, auto_where=None, context=None): """ :arg queue: a :class:`pyopencl.CommandQueue`. - :arg places: a collection or mapping of symbolic names (see also - :class:`pytential.symbolic.execution.GeometryCollection`) to - :class:`meshmode.discretization.Discretization` objects, subclasses - of :class:`pytential.source.PotentialSource` or subclasses of - :class:`pytential.target.TargetBase`. + :arg places: a :class:`pytential.symbolic.execution.GeometryCollection`. + Alternatively, any list or mapping that is a valid argument for its + constructor can also be used. :arg exprs: an array of expressions corresponding to the output block rows of the matrix. May also be a single expression. :arg input_exprs: an array of expressions corresponding to the @@ -597,7 +589,8 @@ def build_matrix(queue, places, exprs, input_exprs, domains=None, # not iterable, wrap in a list input_exprs = [input_exprs] - domains = _prepare_domains(len(input_exprs), places, domains, places.source) + domains = _prepare_domains(len(input_exprs), places, domains, + DEFAULT_SOURCE) from pytential.symbolic.matrix import MatrixBuilder, is_zero nblock_rows = len(exprs) diff --git a/test/test_matrix.py b/test/test_matrix.py index 02eaca0d..3a11bf3a 100644 --- a/test/test_matrix.py +++ b/test/test_matrix.py @@ -209,7 +209,7 @@ def test_p2p_block_builder(ctx_factory, factor, ndim, lpot_id, from pytential.symbolic.execution import _prepare_expr, _prepare_domains places = GeometryCollection(qbx) expr = _prepare_expr(places, op) - domains = _prepare_domains(1, places, None, places.source) + domains = _prepare_domains(1, places, None, DEFAULT_SOURCE) from sumpy.tools import MatrixBlockIndexRanges index_set = MatrixBlockIndexRanges(ctx, tgtindices, srcindices) @@ -376,7 +376,8 @@ def test_build_matrix_where(ctx_factory, where, visualize=False): # build full QBX matrix from pytential.symbolic.execution import build_matrix - mat = build_matrix(queue, qbx, op, u_sym, auto_where=where) + mat = build_matrix(queue, qbx, op, u_sym, + auto_where=where, domains=where[0]) from pytential.symbolic.execution import GeometryCollection places = GeometryCollection(qbx, auto_where=where) -- GitLab From 26d00306de3ba2245014c1b623debb69e51d7463 Mon Sep 17 00:00:00 2001 From: Alex Fikl Date: Wed, 15 Aug 2018 21:35:43 -0500 Subject: [PATCH 29/34] add some copyrights --- pytential/symbolic/execution.py | 5 ++++- pytential/symbolic/matrix.py | 5 ++++- 2 files changed, 8 insertions(+), 2 deletions(-) diff --git a/pytential/symbolic/execution.py b/pytential/symbolic/execution.py index 71300d83..a6528211 100644 --- a/pytential/symbolic/execution.py +++ b/pytential/symbolic/execution.py @@ -1,6 +1,9 @@ from __future__ import division, absolute_import -__copyright__ = "Copyright (C) 2013 Andreas Kloeckner" +__copyright__ = """ +Copyright (C) 2013 Andreas Kloeckner +Copyright (C) 2018 Alexandru Fikl +""" __license__ = """ Permission is hereby granted, free of charge, to any person obtaining a copy diff --git a/pytential/symbolic/matrix.py b/pytential/symbolic/matrix.py index 77cc31ae..9d0c2b2f 100644 --- a/pytential/symbolic/matrix.py +++ b/pytential/symbolic/matrix.py @@ -1,6 +1,9 @@ from __future__ import division, absolute_import -__copyright__ = "Copyright (C) 2015 Andreas Kloeckner" +__copyright__ = """ +Copyright (C) 2015 Andreas Kloeckner +Copyright (C) 2018 Alexandru Fikl +""" __license__ = """ Permission is hereby granted, free of charge, to any person obtaining a copy -- GitLab From fd1eb9ff262718b482813ed7bcabfa1392d5d163 Mon Sep 17 00:00:00 2001 From: Alex Fikl Date: Wed, 22 Aug 2018 18:59:35 -0500 Subject: [PATCH 30/34] execution: store default source and target place ids in GeometryCollection --- pytential/symbolic/execution.py | 9 +++++++-- 1 file changed, 7 insertions(+), 2 deletions(-) diff --git a/pytential/symbolic/execution.py b/pytential/symbolic/execution.py index a6528211..87b32df1 100644 --- a/pytential/symbolic/execution.py +++ b/pytential/symbolic/execution.py @@ -325,7 +325,7 @@ def _prepare_expr(places, expr): from pytential.symbolic.mappers import ( ToTargetTagger, DerivativeBinder) - expr = ToTargetTagger(*places.where)(expr) + expr = ToTargetTagger(*places._default_place_ids)(expr) expr = DerivativeBinder()(expr) for name, place in six.iteritems(places.places): @@ -372,8 +372,13 @@ class GeometryCollection(object): if auto_where is None: source_where, target_where = DEFAULT_SOURCE, DEFAULT_TARGET else: + # NOTE: keeping this here to make sure auto_where unpacks into + # just the two elements source_where, target_where = auto_where - self.where = (source_where, target_where) + + self._default_source_place = source_where + self._default_target_place = target_where + self._auto_place_ids = (source_where, target_where) self.places = {} if isinstance(places, LayerPotentialSourceBase): -- GitLab From 3821f99eb3eb6da99db4441ed95d75755f60c46b Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Andreas=20Kl=C3=B6ckner?= Date: Mon, 27 Aug 2018 17:06:57 -0400 Subject: [PATCH 31/34] Doc/import structure tweaks --- doc/symbolic.rst | 2 ++ pytential/__init__.py | 1 + pytential/symbolic/execution.py | 3 ++- 3 files changed, 5 insertions(+), 1 deletion(-) diff --git a/doc/symbolic.rst b/doc/symbolic.rst index e6de9534..86d580f6 100644 --- a/doc/symbolic.rst +++ b/doc/symbolic.rst @@ -13,6 +13,8 @@ Binding an operator to a discretization .. currentmodule:: pytential +.. autoclass:: GeometryCollection + .. autofunction:: bind PDE operators diff --git a/pytential/__init__.py b/pytential/__init__.py index 1d07499d..8c1f8f2e 100644 --- a/pytential/__init__.py +++ b/pytential/__init__.py @@ -25,6 +25,7 @@ THE SOFTWARE. import numpy as np import pytential.symbolic.primitives as sym +import pytential.symbolic.execution as GeometryCollection from pytential.symbolic.execution import bind from pytools import memoize_on_first_arg diff --git a/pytential/symbolic/execution.py b/pytential/symbolic/execution.py index a6528211..481fc0c8 100644 --- a/pytential/symbolic/execution.py +++ b/pytential/symbolic/execution.py @@ -340,7 +340,8 @@ def _prepare_expr(places, expr): # {{{ bound expression class GeometryCollection(object): - """A collection of geometry-related objects. This class is meant to hold + """A mapping from symbolic identifiers ("place IDs", typically strings) + of layer potential. This class is meant to hold a specific combination of sources and targets and cache any common expressions specific to it, e.g. metric terms, etc. -- GitLab From 1823e1a63bf65cd2593f4f9a3078cd84faa706c3 Mon Sep 17 00:00:00 2001 From: Alex Fikl Date: Mon, 27 Aug 2018 19:30:43 -0500 Subject: [PATCH 32/34] matrix: cleanup tests --- pytential/__init__.py | 2 +- pytential/symbolic/execution.py | 11 +- test/test_matrix.py | 176 +++++++++++++++++++++----------- 3 files changed, 123 insertions(+), 66 deletions(-) diff --git a/pytential/__init__.py b/pytential/__init__.py index 8c1f8f2e..942476f3 100644 --- a/pytential/__init__.py +++ b/pytential/__init__.py @@ -25,7 +25,7 @@ THE SOFTWARE. import numpy as np import pytential.symbolic.primitives as sym -import pytential.symbolic.execution as GeometryCollection +from pytential.symbolic.execution import GeometryCollection from pytential.symbolic.execution import bind from pytools import memoize_on_first_arg diff --git a/pytential/symbolic/execution.py b/pytential/symbolic/execution.py index bd39a1b2..d1e581bb 100644 --- a/pytential/symbolic/execution.py +++ b/pytential/symbolic/execution.py @@ -379,7 +379,7 @@ class GeometryCollection(object): self._default_source_place = source_where self._default_target_place = target_where - self._auto_place_ids = (source_where, target_where) + self._default_place_ids = (source_where, target_where) self.places = {} if isinstance(places, LayerPotentialSourceBase): @@ -407,7 +407,8 @@ class GeometryCollection(object): if not isinstance(lpot, LayerPotentialSourceBase): return lpot - if where is DEFAULT_SOURCE or where is DEFAULT_TARGET: + from pytential.symbolic.primitives import _QBXSource + if not isinstance(where, _QBXSource): where = QBXSourceStage1(where) if isinstance(where, QBXSourceStage1): @@ -575,8 +576,8 @@ def build_matrix(queue, places, exprs, input_exprs, domains=None, :arg domains: a list of discretization identifiers (see 'places') or *None* values indicating the domains on which each component of the solution vector lives. *None* values indicate that the component - is a scalar. If *None*, - :class:`pytential.symbolic.primitives.DEFAULT_SOURCE`, is required + is a scalar. If *None*, *auto_where* or, if it is not provided, + :class:`~pytential.symbolic.primitives.DEFAULT_SOURCE` is required to be a key in :attr:`places`. :arg auto_where: For simple source-to-self or source-to-target evaluations, find 'where' attributes automatically. @@ -599,7 +600,7 @@ def build_matrix(queue, places, exprs, input_exprs, domains=None, input_exprs = [input_exprs] domains = _prepare_domains(len(input_exprs), places, domains, - DEFAULT_SOURCE) + places._default_source_place) from pytential.symbolic.matrix import MatrixBuilder, is_zero nblock_rows = len(exprs) diff --git a/test/test_matrix.py b/test/test_matrix.py index 3a11bf3a..95291380 100644 --- a/test/test_matrix.py +++ b/test/test_matrix.py @@ -31,11 +31,13 @@ import numpy as np import numpy.linalg as la import pyopencl as cl +import pyopencl.array # noqa + from pytools.obj_array import make_obj_array, is_obj_array from sumpy.symbolic import USE_SYMENGINE -from meshmode.mesh.generation import \ - ellipse, NArmedStarfish, make_curve_mesh +from meshmode.mesh.generation import ( # noqa + ellipse, NArmedStarfish, make_curve_mesh, generate_torus) from pytential import bind, sym from pytential.symbolic.primitives import DEFAULT_SOURCE, DEFAULT_TARGET @@ -46,6 +48,77 @@ from pyopencl.tools import ( # noqa pytest_generate_tests_for_pyopencl as pytest_generate_tests) +def _build_qbx_discr(queue, + ndim=2, + nelements=30, + target_order=7, + qbx_order=4, + curve_f=None): + + if curve_f is None: + curve_f = NArmedStarfish(5, 0.25) + + if ndim == 2: + mesh = make_curve_mesh(curve_f, + np.linspace(0, 1, nelements + 1), + target_order) + elif ndim == 3: + mesh = generate_torus(10.0, 2.0, order=target_order) + else: + raise ValueError("unsupported ambient dimension") + + from meshmode.discretization import Discretization + from meshmode.discretization.poly_element import \ + InterpolatoryQuadratureSimplexGroupFactory + from pytential.qbx import QBXLayerPotentialSource + density_discr = Discretization( + queue.context, mesh, + InterpolatoryQuadratureSimplexGroupFactory(target_order)) + + qbx, _ = QBXLayerPotentialSource(density_discr, + fine_order=4 * target_order, + qbx_order=qbx_order, + fmm_order=False).with_refinement() + + return qbx + + +def _build_block_index(discr, nblks=10, factor=1.0): + nnodes = discr.nnodes + max_particles_in_box = nnodes // nblks + + from pytential.linalg.proxy import partition_by_nodes + indices = partition_by_nodes(discr, use_tree=True, + max_nodes_in_box=max_particles_in_box) + + # randomly pick a subset of points + from sumpy.tools import MatrixBlockIndexRanges, BlockIndexRanges + if abs(factor - 1.0) > 1.0e-14: + with cl.CommandQueue(discr.cl_context) as queue: + indices = indices.get(queue) + + indices_ = np.empty(indices.nblocks, dtype=np.object) + for i in range(indices.nblocks): + iidx = indices.block_indices(i) + isize = int(factor * len(iidx)) + isize = max(1, min(isize, len(iidx))) + + indices_[i] = np.sort( + np.random.choice(iidx, size=isize, replace=False)) + + ranges_ = cl.array.to_device(queue, + np.cumsum([0] + [r.shape[0] for r in indices_])) + indices_ = cl.array.to_device(queue, np.hstack(indices_)) + + indices = BlockIndexRanges(discr.cl_context, + indices_.with_queue(None), + ranges_.with_queue(None)) + + indices = MatrixBlockIndexRanges(indices.cl_context, + indices, indices) + + return indices + def _build_op(lpot_id, k=0, @@ -90,6 +163,17 @@ def _build_op(lpot_id, return op, u_sym, knl_kwargs +def _max_block_error(mat, blk, index_set): + error = -np.inf + for i in range(index_set.nblocks): + mat_i = index_set.take(mat, i) + blk_i = index_set.block_take(blk, i) + + error = max(error, la.norm(mat_i - blk_i) / la.norm(mat_i)) + + return error + + @pytest.mark.skipif(USE_SYMENGINE, reason="https://gitlab.tiker.net/inducer/sumpy/issues/25") @pytest.mark.parametrize("k", [0, 42]) @@ -195,15 +279,10 @@ def test_p2p_block_builder(ctx_factory, factor, ndim, lpot_id, from sympy.core.cache import clear_cache clear_cache() - from test_linalg_proxy import _build_qbx_discr, _build_block_index target_order = 2 if ndim == 3 else 7 qbx = _build_qbx_discr(queue, target_order=target_order, ndim=ndim) op, u_sym, _ = _build_op(lpot_id, ndim=ndim) - - srcindices = _build_block_index(qbx.density_discr, - method='nodes', factor=factor) - tgtindices = _build_block_index(qbx.density_discr, - method='nodes', factor=factor) + index_set = _build_block_index(qbx.density_discr, factor=factor) from pytential.symbolic.execution import GeometryCollection from pytential.symbolic.execution import _prepare_expr, _prepare_domains @@ -211,9 +290,6 @@ def test_p2p_block_builder(ctx_factory, factor, ndim, lpot_id, expr = _prepare_expr(places, op) domains = _prepare_domains(1, places, None, DEFAULT_SOURCE) - from sumpy.tools import MatrixBlockIndexRanges - index_set = MatrixBlockIndexRanges(ctx, tgtindices, srcindices) - from pytential.symbolic.matrix import P2PMatrixBuilder mbuilder = P2PMatrixBuilder(queue, dep_expr=u_sym, @@ -257,14 +333,7 @@ def test_p2p_block_builder(ctx_factory, factor, ndim, lpot_id, ax2.set_title('P2PMatrixBuilder') mp.savefig("test_p2p_block_{}d_{:.1f}.png".format(ndim, factor)) - for i in range(index_set.nblocks): - eps = 1.0e-14 * la.norm(index_set.take(mat, i)) - error = la.norm(index_set.block_take(blk, i) - - index_set.take(mat, i)) - - if visualize: - print('block[{:04}]: {:.5e}'.format(i, error)) - assert error < eps + assert _max_block_error(mat, blk, index_set) < 1.0e-14 @pytest.mark.parametrize("factor", [1.0, 0.6]) @@ -279,7 +348,6 @@ def test_qbx_block_builder(ctx_factory, factor, ndim, lpot_id, from sympy.core.cache import clear_cache clear_cache() - from test_linalg_proxy import _build_qbx_discr, _build_block_index target_order = 2 if ndim == 3 else 7 qbx = _build_qbx_discr(queue, target_order=target_order, ndim=ndim) op, u_sym, _ = _build_op(lpot_id, ndim=ndim) @@ -295,11 +363,7 @@ def test_qbx_block_builder(ctx_factory, factor, ndim, lpot_id, density_discr = places.get_discretization(where[0]) from sumpy.tools import MatrixBlockIndexRanges - srcindices = _build_block_index(density_discr, - method='nodes', factor=factor) - tgtindices = _build_block_index(density_discr, - method='nodes', factor=factor) - index_set = MatrixBlockIndexRanges(ctx, tgtindices, srcindices) + index_set = _build_block_index(density_discr, factor=factor) from pytential.symbolic.matrix import NearFieldBlockBuilder mbuilder = NearFieldBlockBuilder(queue, @@ -342,23 +406,16 @@ def test_qbx_block_builder(ctx_factory, factor, ndim, lpot_id, ax2.set_title('NearFieldBlockBuilder') mp.savefig("test_qbx_block_builder.png", dpi=300) - for i in range(index_set.nblocks): - eps = 1.0e-14 * la.norm(index_set.take(mat, i)) - error = la.norm(index_set.block_take(blk, i) - - index_set.take(mat, i)) - - if visualize: - print('block[{:04}]: {:.5e}'.format(i, error)) - assert error < eps + assert _max_block_error(mat, blk, index_set) < 1.0e-14 -@pytest.mark.parametrize('where', +@pytest.mark.parametrize('place_id', [(DEFAULT_SOURCE, DEFAULT_TARGET), (QBXSourceStage1(DEFAULT_SOURCE), QBXSourceStage1(DEFAULT_TARGET)), (QBXSourceQuadStage2(DEFAULT_SOURCE), QBXSourceQuadStage2(DEFAULT_TARGET))]) -def test_build_matrix_where(ctx_factory, where, visualize=False): +def test_build_matrix_places(ctx_factory, place_id, visualize=False): ctx = ctx_factory() queue = cl.CommandQueue(ctx) @@ -366,7 +423,6 @@ def test_build_matrix_where(ctx_factory, where, visualize=False): from sympy.core.cache import clear_cache clear_cache() - from test_linalg_proxy import _build_qbx_discr qbx = _build_qbx_discr(queue, nelements=8, target_order=2, ndim=2, curve_f=partial(ellipse, 1.0)) @@ -374,17 +430,20 @@ def test_build_matrix_where(ctx_factory, where, visualize=False): op, u_sym, _ = _build_op(lpot_id=1, ndim=2, qbx_forced_limit=qbx_forced_limit) + from pytential.symbolic.execution import GeometryCollection + places = GeometryCollection(qbx, auto_where=place_id) + source_discr = places.get_discretization(place_id[0]) + target_discr = places.get_discretization(place_id[1]) + + index_set = _build_block_index(source_discr, factor=0.6) + # build full QBX matrix from pytential.symbolic.execution import build_matrix - mat = build_matrix(queue, qbx, op, u_sym, - auto_where=where, domains=where[0]) - - from pytential.symbolic.execution import GeometryCollection - places = GeometryCollection(qbx, auto_where=where) - source_discr = places.get_discretization(where[0]) - target_discr = places.get_discretization(where[1]) + qbx_mat = build_matrix(queue, qbx, op, u_sym, + auto_where=place_id, domains=place_id[0]) + qbx_mat = qbx_mat.get(queue) - assert mat.shape == (target_discr.nnodes, source_discr.nnodes) + assert qbx_mat.shape == (target_discr.nnodes, source_discr.nnodes) # build full p2p matrix from pytential.symbolic.execution import _prepare_expr @@ -394,43 +453,40 @@ def test_build_matrix_where(ctx_factory, where, visualize=False): mbuilder = P2PMatrixBuilder(queue, dep_expr=u_sym, other_dep_exprs=[], - dep_source=places[where[0]], - dep_discr=places.get_discretization(where[0]), + dep_source=places[place_id[0]], + dep_discr=places.get_discretization(place_id[0]), places=places, context={}) - mat = mbuilder(op) + p2p_mat = mbuilder(op) - assert mat.shape == (target_discr.nnodes, source_discr.nnodes) + assert p2p_mat.shape == (target_discr.nnodes, source_discr.nnodes) # build block qbx and p2p matrices - from test_linalg_proxy import _build_block_index - srcindices = _build_block_index(source_discr, method='nodes', factor=0.6) - tgtindices = _build_block_index(target_discr, method='nodes', factor=0.6) - - from sumpy.tools import MatrixBlockIndexRanges - index_set = MatrixBlockIndexRanges(ctx, tgtindices, srcindices) - from pytential.symbolic.matrix import NearFieldBlockBuilder mbuilder = NearFieldBlockBuilder(queue, dep_expr=u_sym, other_dep_exprs=[], - dep_source=places[where[0]], - dep_discr=places.get_discretization(where[0]), + dep_source=places[place_id[0]], + dep_discr=places.get_discretization(place_id[0]), places=places, index_set=index_set, context={}) mat = mbuilder(op) + if place_id[0] is not DEFAULT_SOURCE: + assert _max_block_error(qbx_mat, mat, index_set.get(queue)) < 1.0e-14 from pytential.symbolic.matrix import FarFieldBlockBuilder mbuilder = FarFieldBlockBuilder(queue, dep_expr=u_sym, other_dep_exprs=[], - dep_source=places[where[0]], - dep_discr=places.get_discretization(where[0]), + dep_source=places[place_id[0]], + dep_discr=places.get_discretization(place_id[0]), places=places, index_set=index_set, - context={}) + context={}, + exclude_self=True) mat = mbuilder(op) + assert _max_block_error(p2p_mat, mat, index_set.get(queue)) < 1.0e-14 if __name__ == "__main__": -- GitLab From b1faa7f885f5545d35f7c5c08649b8edf39c325f Mon Sep 17 00:00:00 2001 From: Alex Fikl Date: Mon, 27 Aug 2018 19:35:31 -0500 Subject: [PATCH 33/34] test_matrix: flake8 fixes --- pytential/__init__.py | 2 +- test/test_matrix.py | 3 +-- 2 files changed, 2 insertions(+), 3 deletions(-) diff --git a/pytential/__init__.py b/pytential/__init__.py index 942476f3..d28e8bdb 100644 --- a/pytential/__init__.py +++ b/pytential/__init__.py @@ -25,7 +25,7 @@ THE SOFTWARE. import numpy as np import pytential.symbolic.primitives as sym -from pytential.symbolic.execution import GeometryCollection +from pytential.symbolic.execution import GeometryCollection # noqa from pytential.symbolic.execution import bind from pytools import memoize_on_first_arg diff --git a/test/test_matrix.py b/test/test_matrix.py index 95291380..d912b2e7 100644 --- a/test/test_matrix.py +++ b/test/test_matrix.py @@ -48,6 +48,7 @@ from pyopencl.tools import ( # noqa pytest_generate_tests_for_pyopencl as pytest_generate_tests) + def _build_qbx_discr(queue, ndim=2, nelements=30, @@ -361,8 +362,6 @@ def test_qbx_block_builder(ctx_factory, factor, ndim, lpot_id, places = GeometryCollection(qbx, auto_where=where) expr = _prepare_expr(places, op) density_discr = places.get_discretization(where[0]) - - from sumpy.tools import MatrixBlockIndexRanges index_set = _build_block_index(density_discr, factor=factor) from pytential.symbolic.matrix import NearFieldBlockBuilder -- GitLab From 4658c114e7f5807b1782b24dc820ae38f8de19e0 Mon Sep 17 00:00:00 2001 From: Andreas Kloeckner Date: Tue, 11 Sep 2018 23:29:47 -0500 Subject: [PATCH 34/34] Doc tweaks --- pytential/source.py | 2 +- pytential/symbolic/execution.py | 11 ++++++++--- 2 files changed, 9 insertions(+), 4 deletions(-) diff --git a/pytential/source.py b/pytential/source.py index e76b97c3..9334dff7 100644 --- a/pytential/source.py +++ b/pytential/source.py @@ -38,7 +38,7 @@ __doc__ = """ class PotentialSource(object): """ - .. method:: preprocess_optemplate(name, discretizations, expr) + .. automethod:: preprocess_optemplate .. method:: op_group_features(expr) diff --git a/pytential/symbolic/execution.py b/pytential/symbolic/execution.py index 78aa9349..976e352f 100644 --- a/pytential/symbolic/execution.py +++ b/pytential/symbolic/execution.py @@ -341,9 +341,13 @@ def _prepare_expr(places, expr): class GeometryCollection(object): """A mapping from symbolic identifiers ("place IDs", typically strings) - of layer potential. This class is meant to hold - a specific combination of sources and targets and cache any common - expressions specific to it, e.g. metric terms, etc. + to 'geometries', where a geometry can be a + :class:`pytential.source.PotentialSource` + or a :class:`pytential.target.TargetBase`. + This class is meant to hold a specific combination of sources and targets + serve to host caches of information derived from them, e.g. FMM trees + of subsets of them, as well as related common subexpressions such as + metric terms. .. method:: __getitem__ .. method:: get_discretization @@ -523,6 +527,7 @@ def bind(places, expr, auto_where=None): if not isinstance(places, GeometryCollection): places = GeometryCollection(places, auto_where=auto_where) + expr = _prepare_expr(places, expr) return BoundExpression(places, expr) -- GitLab