diff --git a/doc/array.rst b/doc/array.rst new file mode 100644 index 0000000000000000000000000000000000000000..a722c4799ace3d3d6f4d10b38333b68404bbbcc0 --- /dev/null +++ b/doc/array.rst @@ -0,0 +1,9 @@ +DOF (Degree-of-Freedom) Arrays +============================== + +.. automodule:: meshmode.dof_array + +Array Contexts +============== + +.. automodule:: meshmode.array_context diff --git a/doc/index.rst b/doc/index.rst index 3a7a2ba71dab133f5f754ef2d213c01f23d97d32..a5e39233edd23cf2ae9fc7a3ec6757d5fffd8fe1 100644 --- a/doc/index.rst +++ b/doc/index.rst @@ -7,6 +7,7 @@ Contents: :maxdepth: 2 mesh + array discretization connection misc diff --git a/examples/plot-connectivity.py b/examples/plot-connectivity.py index 068983ca81672295af65cdf3c82112478edac1a5..1a64312d1d16538f7d5ac62b7e653492c0b94c40 100644 --- a/examples/plot-connectivity.py +++ b/examples/plot-connectivity.py @@ -2,7 +2,8 @@ from __future__ import division import numpy as np # noqa import pyopencl as cl - +from meshmode.array_context import PyOpenCLArrayContext, PyOpenCLArrayContextProfiling +from meshmode.dof_array import thaw order = 4 @@ -10,6 +11,7 @@ order = 4 def main(): cl_ctx = cl.create_some_context() queue = cl.CommandQueue(cl_ctx) + actx = PyOpenCLArrayContextProfiling(queue) from meshmode.mesh.generation import ( # noqa generate_icosphere, generate_icosahedron, @@ -23,13 +25,13 @@ def main(): PolynomialWarpAndBlendGroupFactory discr = Discretization( - cl_ctx, mesh, PolynomialWarpAndBlendGroupFactory(order)) + actx, mesh, PolynomialWarpAndBlendGroupFactory(order)) from meshmode.discretization.visualization import make_visualizer - vis = make_visualizer(queue, discr, order) + vis = make_visualizer(actx, discr, order) vis.write_vtk_file("geometry.vtu", [ - ("f", discr.nodes()[0]), + ("f", thaw(actx, discr.nodes()[0])), ]) from meshmode.discretization.visualization import \ diff --git a/examples/simple-dg.py b/examples/simple-dg.py index 8fe6f70600c3230bca4240b3c0808bf853ea650b..b945687758f172194d859560bc1723064d56971b 100644 --- a/examples/simple-dg.py +++ b/examples/simple-dg.py @@ -27,14 +27,13 @@ import numpy as np import numpy.linalg as la # noqa import pyopencl as cl import pyopencl.array as cla # noqa -import pyopencl.clmath as clmath from pytools import memoize_method, memoize_in from pytools.obj_array import ( - join_fields, make_obj_array, - with_object_array_or_scalar, - is_obj_array) -import loopy as lp + flat_obj_array, make_obj_array, + obj_array_vectorize) from meshmode.mesh import BTAG_ALL, BTAG_NONE # noqa +from meshmode.dof_array import DOFArray, freeze, thaw +from meshmode.array_context import PyOpenCLArrayContext, make_loopy_program # Features lost vs. https://github.com/inducer/grudge: @@ -44,42 +43,39 @@ from meshmode.mesh import BTAG_ALL, BTAG_NONE # noqa # - distributed-memory -def with_queue(queue, ary): - return with_object_array_or_scalar( - lambda x: x.with_queue(queue), ary) - - -def without_queue(ary): - return with_queue(None, ary) - - # {{{ discretization -def parametrization_derivative(queue, discr): +def parametrization_derivative(actx, discr): + thawed_nodes = thaw(actx, discr.nodes()) + result = np.zeros((discr.ambient_dim, discr.dim), dtype=object) for iambient in range(discr.ambient_dim): for idim in range(discr.dim): result[iambient, idim] = discr.num_reference_derivative( - queue, (idim,), discr.nodes()[iambient]) + (idim,), thawed_nodes[iambient]) return result class DGDiscretization: - def __init__(self, cl_ctx, mesh, order): + def __init__(self, actx, mesh, order): self.order = order from meshmode.discretization import Discretization from meshmode.discretization.poly_element import \ PolynomialWarpAndBlendGroupFactory self.group_factory = PolynomialWarpAndBlendGroupFactory(order=order) - self.volume_discr = Discretization(cl_ctx, mesh, self.group_factory) + self.volume_discr = Discretization(actx, mesh, self.group_factory) assert self.volume_discr.dim == 2 @property - def cl_context(self): - return self.volume_discr.cl_context + def _setup_actx(self): + return self.volume_discr._setup_actx + + @property + def array_context(self): + return self.volume_discr.array_context @property def dim(self): @@ -91,6 +87,7 @@ class DGDiscretization: def boundary_connection(self, boundary_tag): from meshmode.discretization.connection import make_face_restriction return make_face_restriction( + self.volume_discr._setup_actx, self.volume_discr, self.group_factory, boundary_tag=boundary_tag) @@ -100,6 +97,7 @@ class DGDiscretization: from meshmode.discretization.connection import ( make_face_restriction, FACE_RESTR_INTERIOR) return make_face_restriction( + self.volume_discr._setup_actx, self.volume_discr, self.group_factory, FACE_RESTR_INTERIOR, @@ -110,13 +108,15 @@ class DGDiscretization: from meshmode.discretization.connection import \ make_opposite_face_connection - return make_opposite_face_connection(self.interior_faces_connection()) + return make_opposite_face_connection( + self._setup_actx, self.interior_faces_connection()) @memoize_method def all_faces_connection(self): from meshmode.discretization.connection import ( make_face_restriction, FACE_RESTR_ALL) return make_face_restriction( + self.volume_discr._setup_actx, self.volume_discr, self.group_factory, FACE_RESTR_ALL, @@ -129,7 +129,7 @@ class DGDiscretization: faces_conn = self.get_connection("vol", where) return make_face_to_all_faces_embedding( - faces_conn, self.get_discr("all_faces")) + self._setup_actx, faces_conn, self.get_discr("all_faces")) def get_connection(self, src, tgt): src_tgt = (src, tgt) @@ -148,11 +148,13 @@ class DGDiscretization: raise ValueError(f"locations '{src}'->'{tgt}' not understood") def interp(self, src, tgt, vec): - if is_obj_array(vec): - return with_object_array_or_scalar( + if (isinstance(vec, np.ndarray) + and vec.dtype.char == "O" + and not isinstance(vec, DOFArray)): + return obj_array_vectorize( lambda el: self.interp(src, tgt, el), vec) - return self.get_connection(src, tgt)(vec.queue, vec) + return self.get_connection(src, tgt)(vec) def get_discr(self, where): if where == "vol": @@ -170,40 +172,35 @@ class DGDiscretization: @memoize_method def parametrization_derivative(self): - with cl.CommandQueue(self.cl_context) as queue: - return without_queue( - parametrization_derivative(queue, self.volume_discr)) + return freeze( + parametrization_derivative(self._setup_actx, self.volume_discr)) @memoize_method def vol_jacobian(self): - with cl.CommandQueue(self.cl_context) as queue: - [a, b], [c, d] = with_queue(queue, self.parametrization_derivative()) - return (a*d-b*c).with_queue(None) + [a, b], [c, d] = thaw(self._setup_actx, self.parametrization_derivative()) + return freeze(a*d-b*c) @memoize_method def inverse_parametrization_derivative(self): - with cl.CommandQueue(self.cl_context) as queue: - [a, b], [c, d] = with_queue(queue, self.parametrization_derivative()) + [a, b], [c, d] = thaw(self._setup_actx, self.parametrization_derivative()) - result = np.zeros((2, 2), dtype=object) - det = a*d-b*c - result[0, 0] = d/det - result[0, 1] = -b/det - result[1, 0] = -c/det - result[1, 1] = a/det + result = np.zeros((2, 2), dtype=object) + det = a*d-b*c + result[0, 0] = d/det + result[0, 1] = -b/det + result[1, 0] = -c/det + result[1, 1] = a/det - return without_queue(result) + return freeze(result) - def zeros(self, queue): - return self.volume_discr.zeros(queue) + def zeros(self, actx): + return self.volume_discr.zeros(actx) def grad(self, vec): ipder = self.inverse_parametrization_derivative() - queue = vec.queue dref = [ - self.volume_discr.num_reference_derivative( - queue, (idim,), vec).with_queue(queue) + self.volume_discr.num_reference_derivative((idim,), vec) for idim in range(self.volume_discr.dim)] return make_obj_array([ @@ -218,22 +215,18 @@ class DGDiscretization: def normal(self, where): bdry_discr = self.get_discr(where) - with cl.CommandQueue(self.cl_context) as queue: - ((a,), (b,)) = with_queue( - queue, parametrization_derivative(queue, bdry_discr)) + ((a,), (b,)) = parametrization_derivative(self._setup_actx, bdry_discr) - nrm = 1/(a**2+b**2)**0.5 - return without_queue(join_fields(b*nrm, -a*nrm)) + nrm = 1/(a**2+b**2)**0.5 + return freeze(flat_obj_array(b*nrm, -a*nrm)) @memoize_method def face_jacobian(self, where): bdry_discr = self.get_discr(where) - with cl.CommandQueue(self.cl_context) as queue: - ((a,), (b,)) = with_queue(queue, - parametrization_derivative(queue, bdry_discr)) + ((a,), (b,)) = parametrization_derivative(self._setup_actx, bdry_discr) - return ((a**2 + b**2)**0.5).with_queue(None) + return freeze((a**2 + b**2)**0.5) @memoize_method def get_inverse_mass_matrix(self, grp, dtype): @@ -242,37 +235,36 @@ class DGDiscretization: grp.basis(), grp.unit_nodes) - with cl.CommandQueue(self.cl_context) as queue: - return (cla.to_device(queue, matrix) - .with_queue(None)) + actx = self._setup_actx + return actx.freeze(actx.from_numpy(matrix)) def inverse_mass(self, vec): - if is_obj_array(vec): - return with_object_array_or_scalar( + if (isinstance(vec, np.ndarray) + and vec.dtype.char == "O" + and not isinstance(vec, DOFArray)): + return obj_array_vectorize( lambda el: self.inverse_mass(el), vec) @memoize_in(self, "elwise_linear_knl") def knl(): - knl = lp.make_kernel( - """{[k,i,j]: - 0<=k all_faces connection -def make_face_to_all_faces_embedding(faces_connection, all_faces_discr, +def make_face_to_all_faces_embedding(actx, faces_connection, all_faces_discr, from_discr=None): """Return a :class:`meshmode.discretization.connection.DiscretizationConnection` @@ -437,61 +432,59 @@ def make_face_to_all_faces_embedding(faces_connection, all_faces_discr, i_faces_grp = 0 - with cl.CommandQueue(vol_discr.cl_context) as queue: - groups = [] - for ivol_grp, vol_grp in enumerate(vol_discr.groups): - batches = [] + groups = [] + for ivol_grp, vol_grp in enumerate(vol_discr.groups): + batches = [] - nfaces = vol_grp.mesh_el_group.nfaces - for iface in range(nfaces): - all_faces_grp = all_faces_discr.groups[i_faces_grp] + nfaces = vol_grp.mesh_el_group.nfaces + for iface in range(nfaces): + all_faces_grp = all_faces_discr.groups[i_faces_grp] - if per_face_groups: - assert len(faces_connection.groups[i_faces_grp].batches) == 1 - else: - assert (len(faces_connection.groups[i_faces_grp].batches) - == nfaces) + if per_face_groups: + assert len(faces_connection.groups[i_faces_grp].batches) == 1 + else: + assert (len(faces_connection.groups[i_faces_grp].batches) + == nfaces) - assert np.array_equal( - from_discr.groups[i_faces_grp].unit_nodes, - all_faces_grp.unit_nodes) + assert np.array_equal( + from_discr.groups[i_faces_grp].unit_nodes, + all_faces_grp.unit_nodes) - # {{{ find src_batch + # {{{ find src_batch - src_batches = faces_connection.groups[i_faces_grp].batches - if per_face_groups: - src_batch, = src_batches - else: - src_batch = src_batches[iface] - del src_batches + src_batches = faces_connection.groups[i_faces_grp].batches + if per_face_groups: + src_batch, = src_batches + else: + src_batch = src_batches[iface] + del src_batches - # }}} + # }}} - if per_face_groups: - to_element_indices = src_batch.from_element_indices - else: - assert all_faces_grp.nelements == nfaces * vol_grp.nelements - - to_element_indices = ( - vol_grp.nelements*iface - + src_batch.from_element_indices.with_queue(queue) - ).with_queue(None) - - batches.append( - InterpolationBatch( - from_group_index=i_faces_grp, - from_element_indices=src_batch.to_element_indices, - to_element_indices=to_element_indices, - result_unit_nodes=all_faces_grp.unit_nodes, - to_element_face=None)) - - is_last_face = iface + 1 == nfaces - if per_face_groups or is_last_face: - groups.append( - DiscretizationConnectionElementGroup(batches=batches)) - batches = [] - - i_faces_grp += 1 + if per_face_groups: + to_element_indices = src_batch.from_element_indices + else: + assert all_faces_grp.nelements == nfaces * vol_grp.nelements + + to_element_indices = actx.freeze( + vol_grp.nelements*iface + + actx.thaw(src_batch.from_element_indices)) + + batches.append( + InterpolationBatch( + from_group_index=i_faces_grp, + from_element_indices=src_batch.to_element_indices, + to_element_indices=to_element_indices, + result_unit_nodes=all_faces_grp.unit_nodes, + to_element_face=None)) + + is_last_face = iface + 1 == nfaces + if per_face_groups or is_last_face: + groups.append( + DiscretizationConnectionElementGroup(batches=batches)) + batches = [] + + i_faces_grp += 1 return DirectDiscretizationConnection( from_discr, diff --git a/meshmode/discretization/connection/opposite_face.py b/meshmode/discretization/connection/opposite_face.py index e084525ba9a6eb5001548f432c25a0a8f861173e..8f83cd4e7e1031a7e319132397899e0cc38ff5cc 100644 --- a/meshmode/discretization/connection/opposite_face.py +++ b/meshmode/discretization/connection/opposite_face.py @@ -26,43 +26,45 @@ from six.moves import range import numpy as np import numpy.linalg as la -import pyopencl as cl -import pyopencl.array # noqa import logging logger = logging.getLogger(__name__) +def freeze_from_numpy(actx, array): + return actx.freeze(actx.from_numpy(array)) + + +def thaw_to_numpy(actx, array): + return actx.to_numpy(actx.thaw(array)) + + # {{{ _make_cross_face_batches -def _make_cross_face_batches(queue, +def _make_cross_face_batches(actx, tgt_bdry_discr, src_bdry_discr, i_tgt_grp, i_src_grp, tgt_bdry_element_indices, src_bdry_element_indices): - def to_dev(ary): - return cl.array.to_device(queue, ary, array_queue=None) from meshmode.discretization.connection.direct import InterpolationBatch if tgt_bdry_discr.dim == 0: yield InterpolationBatch( from_group_index=i_src_grp, - from_element_indices=to_dev(src_bdry_element_indices), - to_element_indices=to_dev(tgt_bdry_element_indices), + from_element_indices=freeze_from_numpy(actx, src_bdry_element_indices), + to_element_indices=freeze_from_numpy(actx, tgt_bdry_element_indices), result_unit_nodes=src_bdry_discr.groups[i_src_grp].unit_nodes, to_element_face=None) return - # FIXME: This should view-then-transfer - # (but PyOpenCL doesn't do non-contiguous transfers for now). - tgt_bdry_nodes = (tgt_bdry_discr.groups[i_tgt_grp] - .view(tgt_bdry_discr.nodes().get(queue=queue)) - [:, tgt_bdry_element_indices]) + tgt_bdry_nodes = np.array([ + thaw_to_numpy(actx, ary[i_tgt_grp])[tgt_bdry_element_indices] + for ary in tgt_bdry_discr.nodes() + ]) - # FIXME: This should view-then-transfer - # (but PyOpenCL doesn't do non-contiguous transfers for now). - src_bdry_nodes = (src_bdry_discr.groups[i_src_grp] - .view(src_bdry_discr.nodes().get(queue=queue)) - [:, src_bdry_element_indices]) + src_bdry_nodes = np.array([ + thaw_to_numpy(actx, ary[i_src_grp])[src_bdry_element_indices] + for ary in src_bdry_discr.nodes() + ]) tol = 1e4 * np.finfo(tgt_bdry_nodes.dtype).eps @@ -255,8 +257,10 @@ def _make_cross_face_batches(queue, from meshmode.discretization.connection.direct import InterpolationBatch yield InterpolationBatch( from_group_index=i_src_grp, - from_element_indices=to_dev(src_bdry_element_indices[close_els]), - to_element_indices=to_dev(tgt_bdry_element_indices[close_els]), + from_element_indices=freeze_from_numpy( + actx, src_bdry_element_indices[close_els]), + to_element_indices=freeze_from_numpy( + actx, tgt_bdry_element_indices[close_els]), result_unit_nodes=template_unit_nodes, to_element_face=None) @@ -276,7 +280,7 @@ def _find_ibatch_for_face(vbc_tgt_grp_batches, iface): return vbc_tgt_grp_face_batch -def _make_bdry_el_lookup_table(queue, connection, igrp): +def _make_bdry_el_lookup_table(actx, connection, igrp): """Given a volume-to-boundary connection as *connection*, return a table of shape ``(from_nelements, nfaces)`` to look up the element number of the boundary element for that face. @@ -289,9 +293,9 @@ def _make_bdry_el_lookup_table(queue, connection, igrp): iel_lookup.fill(-1) for ibatch, batch in enumerate(connection.groups[igrp].batches): - from_element_indices = batch.from_element_indices.get(queue=queue) + from_element_indices = thaw_to_numpy(actx, batch.from_element_indices) iel_lookup[from_element_indices, batch.to_element_face] = \ - batch.to_element_indices.get(queue=queue) + thaw_to_numpy(actx, batch.to_element_indices) return iel_lookup @@ -300,7 +304,7 @@ def _make_bdry_el_lookup_table(queue, connection, igrp): # {{{ make_opposite_face_connection -def make_opposite_face_connection(volume_to_bdry_conn): +def make_opposite_face_connection(actx, volume_to_bdry_conn): """Given a boundary restriction connection *volume_to_bdry_conn*, return a :class:`DirectDiscretizationConnection` that performs data exchange across opposite faces. @@ -323,103 +327,104 @@ def make_opposite_face_connection(volume_to_bdry_conn): # One interpolation batch in this connection corresponds # to a key (i_tgt_grp,) (i_src_grp, i_face_tgt,) - with cl.CommandQueue(vol_discr.cl_context) as queue: - # a list of batches for each group - groups = [[] for i_tgt_grp in range(ngrps)] + # a list of batches for each group + groups = [[] for i_tgt_grp in range(ngrps)] - for i_src_grp in range(ngrps): - src_grp_el_lookup = _make_bdry_el_lookup_table( - queue, volume_to_bdry_conn, i_src_grp) + for i_src_grp in range(ngrps): + src_grp_el_lookup = _make_bdry_el_lookup_table( + actx, volume_to_bdry_conn, i_src_grp) - for i_tgt_grp in range(ngrps): - vbc_tgt_grp_batches = volume_to_bdry_conn.groups[i_tgt_grp].batches + for i_tgt_grp in range(ngrps): + vbc_tgt_grp_batches = volume_to_bdry_conn.groups[i_tgt_grp].batches - adj = vol_mesh.facial_adjacency_groups[i_tgt_grp][i_src_grp] + adj = vol_mesh.facial_adjacency_groups[i_tgt_grp][i_src_grp] - for i_face_tgt in range(vol_mesh.groups[i_tgt_grp].nfaces): - vbc_tgt_grp_face_batch = _find_ibatch_for_face( - vbc_tgt_grp_batches, i_face_tgt) + for i_face_tgt in range(vol_mesh.groups[i_tgt_grp].nfaces): + vbc_tgt_grp_face_batch = _find_ibatch_for_face( + vbc_tgt_grp_batches, i_face_tgt) - # {{{ index wrangling + # {{{ index wrangling - # The elements in the adjacency group will be a subset of - # the elements in the restriction interpolation batch: - # Imagine an inter-group boundary. The volume-to-boundary - # connection will include all faces as targets, whereas - # there will be separate adjacency groups for intra- and - # inter-group connections. + # The elements in the adjacency group will be a subset of + # the elements in the restriction interpolation batch: + # Imagine an inter-group boundary. The volume-to-boundary + # connection will include all faces as targets, whereas + # there will be separate adjacency groups for intra- and + # inter-group connections. - adj_tgt_flags = adj.element_faces == i_face_tgt - adj_els = adj.elements[adj_tgt_flags] - if adj_els.size == 0: - # NOTE: this case can happen for inter-group boundaries - # when all elements are adjacent on the same face - # index, so all other ones will be empty - continue + adj_tgt_flags = adj.element_faces == i_face_tgt + adj_els = adj.elements[adj_tgt_flags] + if adj_els.size == 0: + # NOTE: this case can happen for inter-group boundaries + # when all elements are adjacent on the same face + # index, so all other ones will be empty + continue - vbc_els = vbc_tgt_grp_face_batch.from_element_indices.get(queue) + vbc_els = thaw_to_numpy(actx, + vbc_tgt_grp_face_batch.from_element_indices) - if len(adj_els) == len(vbc_els): - # Same length: assert (below) that the two use the same - # ordering. - vbc_used_els = slice(None) + if len(adj_els) == len(vbc_els): + # Same length: assert (below) that the two use the same + # ordering. + vbc_used_els = slice(None) - else: - # Genuine subset: figure out an index mapping. - vbc_els_sort_idx = np.argsort(vbc_els) - vbc_used_els = vbc_els_sort_idx[np.searchsorted( - vbc_els, adj_els, sorter=vbc_els_sort_idx - )] + else: + # Genuine subset: figure out an index mapping. + vbc_els_sort_idx = np.argsort(vbc_els) + vbc_used_els = vbc_els_sort_idx[np.searchsorted( + vbc_els, adj_els, sorter=vbc_els_sort_idx + )] - assert np.array_equal(vbc_els[vbc_used_els], adj_els) + assert np.array_equal(vbc_els[vbc_used_els], adj_els) - # find to_element_indices + # find to_element_indices - tgt_bdry_element_indices = ( - vbc_tgt_grp_face_batch.to_element_indices - .get(queue=queue)[vbc_used_els]) + tgt_bdry_element_indices = thaw_to_numpy( + actx, + vbc_tgt_grp_face_batch.to_element_indices + )[vbc_used_els] - # find from_element_indices + # find from_element_indices - src_vol_element_indices = adj.neighbors[adj_tgt_flags] - src_element_faces = adj.neighbor_faces[adj_tgt_flags] + src_vol_element_indices = adj.neighbors[adj_tgt_flags] + src_element_faces = adj.neighbor_faces[adj_tgt_flags] - src_bdry_element_indices = src_grp_el_lookup[ - src_vol_element_indices, src_element_faces] + src_bdry_element_indices = src_grp_el_lookup[ + src_vol_element_indices, src_element_faces] - # }}} + # }}} - # {{{ visualization (for debugging) + # {{{ visualization (for debugging) - if 0: - print("TVE", adj.elements[adj_tgt_flags]) - print("TBE", tgt_bdry_element_indices) - print("FVE", src_vol_element_indices) - from meshmode.mesh.visualization import draw_2d_mesh - import matplotlib.pyplot as pt - draw_2d_mesh(vol_discr.mesh, draw_element_numbers=True, - set_bounding_box=True, - draw_vertex_numbers=False, - draw_face_numbers=True, - fill=None) - pt.figure() + if 0: + print("TVE", adj.elements[adj_tgt_flags]) + print("TBE", tgt_bdry_element_indices) + print("FVE", src_vol_element_indices) + from meshmode.mesh.visualization import draw_2d_mesh + import matplotlib.pyplot as pt + draw_2d_mesh(vol_discr.mesh, draw_element_numbers=True, + set_bounding_box=True, + draw_vertex_numbers=False, + draw_face_numbers=True, + fill=None) + pt.figure() - draw_2d_mesh(bdry_discr.mesh, draw_element_numbers=True, - set_bounding_box=True, - draw_vertex_numbers=False, - draw_face_numbers=True, - fill=None) + draw_2d_mesh(bdry_discr.mesh, draw_element_numbers=True, + set_bounding_box=True, + draw_vertex_numbers=False, + draw_face_numbers=True, + fill=None) - pt.show() + pt.show() - # }}} + # }}} - batches = _make_cross_face_batches(queue, - bdry_discr, bdry_discr, - i_tgt_grp, i_src_grp, - tgt_bdry_element_indices, - src_bdry_element_indices) - groups[i_tgt_grp].extend(batches) + batches = _make_cross_face_batches(actx, + bdry_discr, bdry_discr, + i_tgt_grp, i_src_grp, + tgt_bdry_element_indices, + src_bdry_element_indices) + groups[i_tgt_grp].extend(batches) from meshmode.discretization.connection import ( DirectDiscretizationConnection, DiscretizationConnectionElementGroup) @@ -436,7 +441,7 @@ def make_opposite_face_connection(volume_to_bdry_conn): # {{{ make_partition_connection -def make_partition_connection(local_bdry_conn, i_local_part, +def make_partition_connection(actx, local_bdry_conn, i_local_part, remote_bdry, remote_adj_groups, remote_from_elem_faces, remote_from_elem_indices): """ @@ -472,52 +477,50 @@ def make_partition_connection(local_bdry_conn, i_local_part, part_batches = [[] for _ in local_groups] - with cl.CommandQueue(local_bdry_conn.cl_context) as queue: - - for i_remote_grp, adj in enumerate(remote_adj_groups): - indices = (i_local_part == adj.neighbor_partitions) - if not np.any(indices): - # Skip because i_remote_grp is not connected to i_local_part. - continue - i_remote_faces = adj.element_faces[indices] - i_local_meshwide_elems = adj.global_neighbors[indices] - i_local_faces = adj.neighbor_faces[indices] - - i_local_grps = find_group_indices(local_groups, i_local_meshwide_elems) - - for i_local_grp in np.unique(i_local_grps): - - elem_base = local_groups[i_local_grp].element_nr_base - local_el_lookup = _make_bdry_el_lookup_table(queue, - local_bdry_conn, - i_local_grp) - - for i_remote_face in i_remote_faces: - - index_flags = np.logical_and(i_local_grps == i_local_grp, - i_remote_faces == i_remote_face) - if not np.any(index_flags): - continue - - remote_bdry_indices = None - for idxs, face in zip(remote_from_elem_indices[i_remote_grp], - remote_from_elem_faces[i_remote_grp]): - if face == i_remote_face: - remote_bdry_indices = idxs - break - assert remote_bdry_indices is not None - - elems = i_local_meshwide_elems[index_flags] - elem_base - faces = i_local_faces[index_flags] - local_bdry_indices = local_el_lookup[elems, faces] - - batches = _make_cross_face_batches(queue, - local_bdry, remote_bdry, - i_local_grp, i_remote_grp, - local_bdry_indices, - remote_bdry_indices) - - part_batches[i_local_grp].extend(batches) + for i_remote_grp, adj in enumerate(remote_adj_groups): + indices = (i_local_part == adj.neighbor_partitions) + if not np.any(indices): + # Skip because i_remote_grp is not connected to i_local_part. + continue + i_remote_faces = adj.element_faces[indices] + i_local_meshwide_elems = adj.global_neighbors[indices] + i_local_faces = adj.neighbor_faces[indices] + + i_local_grps = find_group_indices(local_groups, i_local_meshwide_elems) + + for i_local_grp in np.unique(i_local_grps): + + elem_base = local_groups[i_local_grp].element_nr_base + local_el_lookup = _make_bdry_el_lookup_table(actx, + local_bdry_conn, + i_local_grp) + + for i_remote_face in i_remote_faces: + + index_flags = np.logical_and(i_local_grps == i_local_grp, + i_remote_faces == i_remote_face) + if not np.any(index_flags): + continue + + remote_bdry_indices = None + for idxs, face in zip(remote_from_elem_indices[i_remote_grp], + remote_from_elem_faces[i_remote_grp]): + if face == i_remote_face: + remote_bdry_indices = idxs + break + assert remote_bdry_indices is not None + + elems = i_local_meshwide_elems[index_flags] - elem_base + faces = i_local_faces[index_flags] + local_bdry_indices = local_el_lookup[elems, faces] + + batches = _make_cross_face_batches(actx, + local_bdry, remote_bdry, + i_local_grp, i_remote_grp, + local_bdry_indices, + remote_bdry_indices) + + part_batches[i_local_grp].extend(batches) return DirectDiscretizationConnection( from_discr=remote_bdry, diff --git a/meshmode/discretization/connection/projection.py b/meshmode/discretization/connection/projection.py index fae3f9b7435b526cd7342090bd31012885945ff9..575b07aa5a945ee346d7f67af03da001455a029a 100644 --- a/meshmode/discretization/connection/projection.py +++ b/meshmode/discretization/connection/projection.py @@ -23,12 +23,13 @@ THE SOFTWARE. """ import numpy as np -import pyopencl as cl -import pyopencl.array # noqa -from loopy.version import MOST_RECENT_LANGUAGE_VERSION -from pytools import memoize_method, memoize_in +from pytools import keyed_memoize_method, memoize_in +import loopy as lp + +from meshmode.array_context import make_loopy_program +from meshmode.dof_array import DOFArray from meshmode.discretization.connection.direct import ( DiscretizationConnection, DirectDiscretizationConnection) @@ -78,8 +79,8 @@ class L2ProjectionInverseDiscretizationConnection(DiscretizationConnection): to_discr=self.conn.from_discr, is_surjective=is_surjective) - @memoize_method - def _batch_weights(self): + @keyed_memoize_method(key=lambda actx: ()) + def _batch_weights(self, actx): """Computes scaled quadrature weights for each interpolation batch in :attr:`conn`. The quadrature weights can be used to integrate over child elements in the domain of the parent element, by a change of @@ -111,26 +112,32 @@ class L2ProjectionInverseDiscretizationConnection(DiscretizationConnection): mat = grp.diff_matrices()[iaxis] jac[iaxis] = mat.dot(batch.result_unit_nodes.T) - weights[igrp, ibatch] = det(jac) * grp.weights + weights[igrp, ibatch] = actx.freeze(actx.from_numpy( + det(jac) * grp.weights)) return weights - def __call__(self, queue, vec): - @memoize_in(self, "conn_projection_knl") + def __call__(self, vec): + if not isinstance(vec, DOFArray): + raise TypeError("non-array passed to discretization connection") + + actx = vec.array_context + + @memoize_in(actx, (L2ProjectionInverseDiscretizationConnection, + "conn_projection_knl")) def kproj(): - import loopy as lp - knl = lp.make_kernel([ - "{[k]: 0 <= k < nelements}", - "{[j]: 0 <= j < n_from_nodes}" + return make_loopy_program([ + "{[iel]: 0 <= iel < nelements}", + "{[idof_quad]: 0 <= idof_quad < n_from_nodes}" ], """ - for k - <> element_dot = \ - sum(j, vec[from_element_indices[k], j] * \ - basis[j] * weights[j]) + for iel + <> element_dot = sum(idof_quad, + vec[from_element_indices[iel], idof_quad] + * basis[idof_quad] * weights[idof_quad]) - result[to_element_indices[k], ibasis] = \ - result[to_element_indices[k], ibasis] + element_dot + result[to_element_indices[iel], ibasis] = \ + result[to_element_indices[iel], ibasis] + element_dot end """, [ @@ -148,21 +155,18 @@ class L2ProjectionInverseDiscretizationConnection(DiscretizationConnection): lp.ValueArg("ibasis", np.int32), '...' ], - name="conn_projection_knl", - lang_version=MOST_RECENT_LANGUAGE_VERSION) - - return knl + name="conn_projection_knl") - @memoize_in(self, "conn_evaluation_knl") + @memoize_in(actx, (L2ProjectionInverseDiscretizationConnection, + "conn_evaluation_knl")) def keval(): - import loopy as lp - knl = lp.make_kernel([ - "{[k]: 0 <= k < nelements}", - "{[j]: 0 <= j < n_to_nodes}" + return make_loopy_program([ + "{[iel]: 0 <= iel < nelements}", + "{[idof]: 0 <= idof < n_to_nodes}" ], """ - result[k, j] = result[k, j] + \ - coefficients[k, ibasis] * basis[j] + result[iel, idof] = result[iel, idof] + \ + coefficients[iel, ibasis] * basis[idof] """, [ lp.GlobalArg("coefficients", None, @@ -170,55 +174,49 @@ class L2ProjectionInverseDiscretizationConnection(DiscretizationConnection): lp.ValueArg("ibasis", np.int32), '...' ], - name="conn_evaluate_knl", - default_offset=lp.auto, - lang_version=MOST_RECENT_LANGUAGE_VERSION) - - return knl - - if not isinstance(vec, cl.array.Array): - raise TypeError("non-array passed to discretization connection") - - if vec.shape != (self.from_discr.nnodes,): - raise ValueError("invalid shape of incoming resampling data") + name="conn_evaluate_knl") # compute weights on each refinement of the reference element - weights = self._batch_weights() + weights = self._batch_weights(actx) # perform dot product (on reference element) to get basis coefficients - c = self.to_discr.zeros(queue, dtype=vec.dtype) + c = self.to_discr.zeros(actx, dtype=vec.entry_dtype) + for igrp, (tgrp, cgrp) in enumerate( zip(self.to_discr.groups, self.conn.groups)): for ibatch, batch in enumerate(cgrp.batches): sgrp = self.from_discr.groups[batch.from_group_index] for ibasis, basis_fn in enumerate(sgrp.basis()): - basis = basis_fn(batch.result_unit_nodes).flatten() + basis = actx.from_numpy( + basis_fn(batch.result_unit_nodes).flatten()) # NOTE: batch.*_element_indices are reversed here because # they are from the original forward connection, but # we are going in reverse here. a bit confusing, but # saves on recreating the connection groups and batches. - kproj()(queue, + actx.call_loopy(kproj(), ibasis=ibasis, - vec=sgrp.view(vec), + vec=vec[sgrp.index], basis=basis, weights=weights[igrp, ibatch], - result=tgrp.view(c), + result=c[igrp], from_element_indices=batch.to_element_indices, to_element_indices=batch.from_element_indices) # evaluate at unit_nodes to get the vector on to_discr - result = self.to_discr.zeros(queue, dtype=vec.dtype) + result = self.to_discr.zeros(actx, dtype=vec.entry_dtype) for igrp, grp in enumerate(self.to_discr.groups): for ibasis, basis_fn in enumerate(grp.basis()): - basis = basis_fn(grp.unit_nodes).flatten() + basis = actx.from_numpy( + basis_fn(grp.unit_nodes).flatten()) - keval()(queue, + actx.call_loopy( + keval(), ibasis=ibasis, - result=grp.view(result), + result=result[grp.index], basis=basis, - coefficients=grp.view(c)) + coefficients=c[grp.index]) return result diff --git a/meshmode/discretization/connection/refinement.py b/meshmode/discretization/connection/refinement.py index 9a2d60d890f6b436e72151d66171d819cdce6a98..55026e59a453c5a6fd1738c381ac313863242d15 100644 --- a/meshmode/discretization/connection/refinement.py +++ b/meshmode/discretization/connection/refinement.py @@ -24,8 +24,6 @@ THE SOFTWARE. """ import numpy as np -import pyopencl as cl -import pyopencl.array # noqa import logging logger = logging.getLogger(__name__) @@ -36,7 +34,7 @@ from pytools import log_process # {{{ Build interpolation batches for group def _build_interpolation_batches_for_group( - queue, group_idx, coarse_discr_group, fine_discr_group, record): + actx, group_idx, coarse_discr_group, fine_discr_group, record): r""" To map between discretizations, we sort each of the fine mesh elements into an interpolation batch. Which batch they go @@ -104,8 +102,8 @@ def _build_interpolation_batches_for_group( continue yield InterpolationBatch( from_group_index=group_idx, - from_element_indices=cl.array.to_device(queue, np.asarray(from_bin)), - to_element_indices=cl.array.to_device(queue, np.asarray(to_bin)), + from_element_indices=actx.from_numpy(np.asarray(from_bin)), + to_element_indices=actx.from_numpy(np.asarray(to_bin)), result_unit_nodes=unit_nodes, to_element_face=None) @@ -113,7 +111,7 @@ def _build_interpolation_batches_for_group( @log_process(logger) -def make_refinement_connection(refiner, coarse_discr, group_factory): +def make_refinement_connection(actx, refiner, coarse_discr, group_factory): """Return a :class:`meshmode.discretization.connection.DiscretizationConnection` connecting `coarse_discr` to a discretization on the fine mesh. @@ -142,21 +140,20 @@ def make_refinement_connection(refiner, coarse_discr, group_factory): from meshmode.discretization import Discretization fine_discr = Discretization( - coarse_discr.cl_context, + actx, fine_mesh, group_factory, real_dtype=coarse_discr.real_dtype) groups = [] - with cl.CommandQueue(fine_discr.cl_context) as queue: - for group_idx, (coarse_discr_group, fine_discr_group, record) in \ - enumerate(zip(coarse_discr.groups, fine_discr.groups, - refiner.group_refinement_records)): - groups.append( - DiscretizationConnectionElementGroup( - list(_build_interpolation_batches_for_group( - queue, group_idx, coarse_discr_group, - fine_discr_group, record)))) + for group_idx, (coarse_discr_group, fine_discr_group, record) in \ + enumerate(zip(coarse_discr.groups, fine_discr.groups, + refiner.group_refinement_records)): + groups.append( + DiscretizationConnectionElementGroup( + list(_build_interpolation_batches_for_group( + actx, group_idx, coarse_discr_group, + fine_discr_group, record)))) return DirectDiscretizationConnection( from_discr=coarse_discr, diff --git a/meshmode/discretization/connection/same_mesh.py b/meshmode/discretization/connection/same_mesh.py index e7b781cd8e2908e1ef9b5af7da202a009931a256..355d534944e905974ca61ec210f64bb48e91bfcd 100644 --- a/meshmode/discretization/connection/same_mesh.py +++ b/meshmode/discretization/connection/same_mesh.py @@ -23,13 +23,11 @@ THE SOFTWARE. """ import numpy as np -import pyopencl as cl -import pyopencl.array # noqa # {{{ same-mesh constructor -def make_same_mesh_connection(to_discr, from_discr): +def make_same_mesh_connection(actx, to_discr, from_discr): from meshmode.discretization.connection.direct import ( InterpolationBatch, DiscretizationConnectionElementGroup, @@ -39,23 +37,22 @@ def make_same_mesh_connection(to_discr, from_discr): raise ValueError("from_discr and to_discr must be based on " "the same mesh") - assert to_discr.cl_context == from_discr.cl_context - - with cl.CommandQueue(to_discr.cl_context) as queue: - groups = [] - for igrp, (fgrp, tgrp) in enumerate(zip(from_discr.groups, to_discr.groups)): - all_elements = cl.array.arange(queue, - fgrp.nelements, - dtype=np.intp).with_queue(None) - ibatch = InterpolationBatch( - from_group_index=igrp, - from_element_indices=all_elements, - to_element_indices=all_elements, - result_unit_nodes=tgrp.unit_nodes, - to_element_face=None) - - groups.append( - DiscretizationConnectionElementGroup([ibatch])) + groups = [] + for igrp, (fgrp, tgrp) in enumerate(zip(from_discr.groups, to_discr.groups)): + all_elements = actx.freeze( + actx.from_numpy( + np.arange( + fgrp.nelements, + dtype=np.intp))) + ibatch = InterpolationBatch( + from_group_index=igrp, + from_element_indices=all_elements, + to_element_indices=all_elements, + result_unit_nodes=tgrp.unit_nodes, + to_element_face=None) + + groups.append( + DiscretizationConnectionElementGroup([ibatch])) return DirectDiscretizationConnection( from_discr, to_discr, groups, diff --git a/meshmode/discretization/poly_element.py b/meshmode/discretization/poly_element.py index 436b1fb1229e2dc11bc56faa280cf0185763735c..91fccdf645862b13233df5e639b4ed02c29cf53d 100644 --- a/meshmode/discretization/poly_element.py +++ b/meshmode/discretization/poly_element.py @@ -53,15 +53,6 @@ Group factories .. autoclass:: LegendreGaussLobattoTensorProductGroupFactory """ -# FIXME Most of the loopy kernels will break as soon as we start using multiple -# element groups. That's because then the dimension-to-dimension stride will no -# longer just be the long axis of the array, but something entirely -# independent. The machinery for this on the loopy end is there, in the form -# of the "stride:auto" dim tag, it just needs to be pushed through all the -# kernels. Fortunately, this will fail in an obvious and noisy way, because -# loopy sees strides that it doesn't expect and complains. - - from meshmode.discretization import ElementGroupBase, InterpolatoryElementGroupBase @@ -300,12 +291,12 @@ class HomogeneousOrderBasedGroupFactory(ElementGroupFactory): def __init__(self, order): self.order = order - def __call__(self, mesh_el_group, node_nr_base): + def __call__(self, mesh_el_group, index): if not isinstance(mesh_el_group, self.mesh_group_class): raise TypeError("only mesh element groups of type '%s' " "are supported" % self.mesh_group_class.__name__) - return self.group_class(mesh_el_group, self.order, node_nr_base) + return self.group_class(mesh_el_group, self.order, index) class OrderAndTypeBasedGroupFactory(ElementGroupFactory): @@ -314,7 +305,7 @@ class OrderAndTypeBasedGroupFactory(ElementGroupFactory): self.simplex_group_class = simplex_group_class self.tensor_product_group_class = tensor_product_group_class - def __call__(self, mesh_el_group, node_nr_base): + def __call__(self, mesh_el_group, index): if isinstance(mesh_el_group, _MeshSimplexElementGroup): group_class = self.simplex_group_class elif isinstance(mesh_el_group, _MeshTensorProductElementGroup): @@ -325,7 +316,7 @@ class OrderAndTypeBasedGroupFactory(ElementGroupFactory): _MeshSimplexElementGroup.__name__, _MeshTensorProductElementGroup.__name__)) - return group_class(mesh_el_group, self.order, node_nr_base) + return group_class(mesh_el_group, self.order, index) # {{{ group factories for simplices diff --git a/meshmode/discretization/visualization.py b/meshmode/discretization/visualization.py index e985f3ffcd9f077a2c6dcb97b6e94e4c8147c49c..24fe92f675293ef6ab048553778ac50a394d4173 100644 --- a/meshmode/discretization/visualization.py +++ b/meshmode/discretization/visualization.py @@ -1,5 +1,3 @@ -from __future__ import division, absolute_import - __copyright__ = "Copyright (C) 2014 Andreas Kloeckner" __license__ = """ @@ -22,10 +20,10 @@ OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN THE SOFTWARE. """ -from six.moves import range import numpy as np from pytools import memoize_method, Record -import pyopencl as cl +from meshmode.dof_array import DOFArray, flatten, thaw + __doc__ = """ @@ -40,29 +38,28 @@ __doc__ = """ # {{{ visualizer def separate_by_real_and_imag(data, real_only): - for name, field in data: - from pytools.obj_array import log_shape, is_obj_array - ls = log_shape(field) + # This function is called on numpy data that has already been + # merged into a single vector. - if is_obj_array(field): - assert len(ls) == 1 + for name, field in data: + if isinstance(field, np.ndarray) and field.dtype.char == "O": + assert len(field.shape) == 1 from pytools.obj_array import ( - oarray_real_copy, oarray_imag_copy, - with_object_array_or_scalar) + obj_array_real_copy, obj_array_imag_copy, + obj_array_vectorize) if field[0].dtype.kind == "c": if real_only: yield (name, - with_object_array_or_scalar(oarray_real_copy, field)) + obj_array_vectorize(obj_array_real_copy, field)) else: yield (name+"_r", - with_object_array_or_scalar(oarray_real_copy, field)) + obj_array_vectorize(obj_array_real_copy, field)) yield (name+"_i", - with_object_array_or_scalar(oarray_imag_copy, field)) + obj_array_vectorize(obj_array_imag_copy, field)) else: yield (name, field) else: - # ls == () if field.dtype.kind == "c": yield (name+"_r", field.real.copy()) yield (name+"_i", field.imag.copy()) @@ -115,17 +112,29 @@ class Visualizer(object): self.element_shrink_factor = element_shrink_factor - def _resample_and_get(self, queue, vec): - from pytools.obj_array import with_object_array_or_scalar - - def resample_and_get_one(fld): - from numbers import Number - if isinstance(fld, Number): - return np.ones(self.connection.to_discr.nnodes) * fld - else: - return self.connection(queue, fld).get(queue=queue) + def _resample_to_numpy(self, vec): + if (isinstance(vec, np.ndarray) + and vec.dtype.char == "O" + and not isinstance(vec, DOFArray)): + from pytools.obj_array import obj_array_vectorize + return obj_array_vectorize(self._resample_to_numpy, vec) + + from numbers import Number + if isinstance(vec, Number): + nnodes = sum(grp.ndofs for grp in self.connection.to_discr.groups) + return np.ones(nnodes) * vec + else: + resampled = self.connection(vec) + actx = resampled.array_context + return actx.to_numpy(flatten(resampled)) - return with_object_array_or_scalar(resample_and_get_one, vec) + @memoize_method + def _vis_nodes(self): + actx = self.vis_discr._setup_actx + return np.array([ + actx.to_numpy(flatten(thaw(actx, ary))) + for ary in self.vis_discr.nodes() + ]) # {{{ vis sub-element connectivity @@ -148,14 +157,15 @@ class Visualizer(object): VTK_QUAD, VTK_HEXAHEDRON) subel_nr_base = 0 + node_nr_base = 0 for group in self.vis_discr.groups: if isinstance(group.mesh_el_group, SimplexElementGroup): node_tuples = list(gnitstam(group.order, group.dim)) - from modepy.tools import submesh + from modepy.tools import simplex_submesh el_connectivity = np.array( - submesh(node_tuples), + simplex_submesh(node_tuples), dtype=np.intp) vtk_cell_type = { 1: VTK_LINE, @@ -202,10 +212,10 @@ class Visualizer(object): raise NotImplementedError("visualization for element groups " "of type '%s'" % type(group.mesh_el_group).__name__) - assert len(node_tuples) == group.nunit_nodes + assert len(node_tuples) == group.nunit_dofs vis_connectivity = ( - group.node_nr_base + np.arange( - 0, group.nelements*group.nunit_nodes, group.nunit_nodes + node_nr_base + np.arange( + 0, group.nelements*group.nunit_dofs, group.nunit_dofs )[:, np.newaxis, np.newaxis] + el_connectivity).astype(np.intp) @@ -216,6 +226,7 @@ class Visualizer(object): result.append(vgrp) subel_nr_base += vgrp.nsubelements + node_nr_base += group.ndofs return result @@ -228,10 +239,8 @@ class Visualizer(object): do_show = kwargs.pop("do_show", True) - with cl.CommandQueue(self.vis_discr.cl_context) as queue: - nodes = self.vis_discr.nodes().with_queue(queue).get() - - field = self._resample_and_get(queue, field) + nodes = self._vis_nodes() + field = self._resample_to_numpy(field) assert nodes.shape[0] == self.vis_discr.ambient_dim #mlab.points3d(nodes[0], nodes[1], 0*nodes[0]) @@ -286,12 +295,10 @@ class Visualizer(object): AppendedDataXMLGenerator, VF_LIST_OF_COMPONENTS) - with cl.CommandQueue(self.vis_discr.cl_context) as queue: - nodes = self.vis_discr.nodes().with_queue(queue).get() - - names_and_fields = [ - (name, self._resample_and_get(queue, fld)) - for name, fld in names_and_fields] + nodes = self._vis_nodes() + names_and_fields = [ + (name, self._resample_to_numpy(fld)) + for name, fld in names_and_fields] vc_groups = self._vis_connectivity() @@ -318,8 +325,12 @@ class Visualizer(object): + (1-self.element_shrink_factor) * el_centers[:, :, np.newaxis]) + if len(self.vis_discr.groups) != 1: + raise NotImplementedError("visualization with multiple " + "element groups") + grid = UnstructuredGrid( - (self.vis_discr.nnodes, + (sum(grp.ndofs for grp in self.vis_discr.groups), DataArray("points", nodes.reshape(self.vis_discr.ambient_dim, -1), vector_format=VF_LIST_OF_COMPONENTS)), @@ -362,10 +373,8 @@ class Visualizer(object): vmax = kwargs.pop("vmax", None) norm = kwargs.pop("norm", None) - with cl.CommandQueue(self.vis_discr.cl_context) as queue: - nodes = self.vis_discr.nodes().with_queue(queue).get() - - field = self._resample_and_get(queue, field) + nodes = self._vis_nodes() + field = self._resample_to_numpy(field) assert nodes.shape[0] == self.vis_discr.ambient_dim @@ -420,14 +429,14 @@ class Visualizer(object): # }}} -def make_visualizer(queue, discr, vis_order, element_shrink_factor=None): +def make_visualizer(actx, discr, vis_order, element_shrink_factor=None): from meshmode.discretization import Discretization from meshmode.discretization.poly_element import ( PolynomialWarpAndBlendElementGroup, LegendreGaussLobattoTensorProductElementGroup, OrderAndTypeBasedGroupFactory) vis_discr = Discretization( - discr.cl_context, discr.mesh, + actx, discr.mesh, OrderAndTypeBasedGroupFactory( vis_order, simplex_group_class=PolynomialWarpAndBlendElementGroup, @@ -438,7 +447,7 @@ def make_visualizer(queue, discr, vis_order, element_shrink_factor=None): make_same_mesh_connection return Visualizer( - make_same_mesh_connection(vis_discr, discr), + make_same_mesh_connection(actx, vis_discr, discr), element_shrink_factor=element_shrink_factor) # }}} @@ -453,16 +462,18 @@ def draw_curve(discr): plt.plot(mesh.vertices[0], mesh.vertices[1], "o") color = plt.cm.rainbow(np.linspace(0, 1, len(discr.groups))) - with cl.CommandQueue(discr.cl_context) as queue: - for i, group in enumerate(discr.groups): - group_nodes = group.view(discr.nodes()).get(queue=queue) - artist_handles = plt.plot( - group_nodes[0].T, - group_nodes[1].T, "-x", - color=color[i]) - - if artist_handles: - artist_handles[0].set_label("Group %d" % i) + for igrp, group in enumerate(discr.groups): + group_nodes = np.array([ + discr._setup_actx.to_numpy(discr.nodes()[iaxis][igrp]) + for iaxis in range(discr.ambient_dim) + ]) + artist_handles = plt.plot( + group_nodes[0].T, + group_nodes[1].T, "-x", + color=color[igrp]) + + if artist_handles: + artist_handles[0].set_label("Group %d" % igrp) # }}} diff --git a/meshmode/distributed.py b/meshmode/distributed.py index 54135fb061097ed8064b696b1b4b778e70dd02c4..bbe24bb0511b072a89192f2703c2e232f7dfcf7c 100644 --- a/meshmode/distributed.py +++ b/meshmode/distributed.py @@ -128,15 +128,14 @@ class MPIBoundaryCommSetupHelper(object): .. automethod:: __call__ .. automethod:: is_ready """ - def __init__(self, mpi_comm, queue, local_bdry_conn, i_remote_part, + def __init__(self, mpi_comm, actx, local_bdry_conn, i_remote_part, bdry_grp_factory): """ :arg mpi_comm: A :class:`MPI.Intracomm` - :arg queue: :arg i_remote_part: The part number of the remote partition """ self.mpi_comm = mpi_comm - self.queue = queue + self.array_context = actx self.i_local_part = mpi_comm.Get_rank() self.i_remote_part = i_remote_part self.local_bdry_conn = local_bdry_conn @@ -151,9 +150,11 @@ class MPIBoundaryCommSetupHelper(object): for i in range(len(local_mesh.groups))] local_to_elem_faces = [[batch.to_element_face for batch in grp_batches] for grp_batches in local_batches] - local_to_elem_indices = [[batch.to_element_indices.get(queue=self.queue) - for batch in grp_batches] - for grp_batches in local_batches] + local_to_elem_indices = [ + [ + self.array_context.to_numpy(batch.to_element_indices) + for batch in grp_batches] + for grp_batches in local_batches] local_data = {'bdry_mesh': local_bdry.mesh, 'adj': local_adj_groups, @@ -190,7 +191,7 @@ class MPIBoundaryCommSetupHelper(object): from meshmode.discretization import Discretization remote_bdry_mesh = remote_data['bdry_mesh'] - remote_bdry = Discretization(self.queue.context, remote_bdry_mesh, + remote_bdry = Discretization(self.array_context, remote_bdry_mesh, self.bdry_grp_factory) remote_adj_groups = remote_data['adj'] remote_to_elem_faces = remote_data['to_elem_faces'] @@ -198,12 +199,11 @@ class MPIBoundaryCommSetupHelper(object): # Connect local_mesh to remote_mesh from meshmode.discretization.connection import make_partition_connection - remote_to_local_bdry_conn = make_partition_connection(self.local_bdry_conn, - self.i_local_part, - remote_bdry, - remote_adj_groups, - remote_to_elem_faces, - remote_to_elem_indices) + remote_to_local_bdry_conn = make_partition_connection( + self.array_context, self.local_bdry_conn, self.i_local_part, + remote_bdry, remote_adj_groups, remote_to_elem_faces, + remote_to_elem_indices) + self.send_req.wait() return remote_to_local_bdry_conn diff --git a/meshmode/dof_array.py b/meshmode/dof_array.py new file mode 100644 index 0000000000000000000000000000000000000000..d04dedc49925c780fec300a023342fd9454c1e16 --- /dev/null +++ b/meshmode/dof_array.py @@ -0,0 +1,244 @@ +__copyright__ = "Copyright (C) 2020 Andreas Kloeckner" + +__license__ = """ +Permission is hereby granted, free of charge, to any person obtaining a copy +of this software and associated documentation files (the "Software"), to deal +in the Software without restriction, including without limitation the rights +to use, copy, modify, merge, publish, distribute, sublicense, and/or sell +copies of the Software, and to permit persons to whom the Software is +furnished to do so, subject to the following conditions: + +The above copyright notice and this permission notice shall be included in +all copies or substantial portions of the Software. + +THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR +IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY, +FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL THE +AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER +LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM, +OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN +THE SOFTWARE. +""" + +import numpy as np +from typing import Optional, TYPE_CHECKING +from functools import partial + +from pytools import single_valued, memoize_in +from pytools.obj_array import obj_array_vectorize + +from meshmode.array_context import ArrayContext, make_loopy_program + +if TYPE_CHECKING: + from meshmode.discretization import Discretization as _Discretization + + +__doc__ = """ +.. autoclass:: DOFArray +.. autofunction:: thaw +.. autofunction:: freeze +.. autofunction:: flatten +.. autofunction:: unflatten +.. autofunction:: flat_norm +""" + + +# {{{ DOFArray + +class DOFArray(np.ndarray): + """This array type is a subclass of :class:`numpy.ndarray` intended to hold + degree-of-freedom arrays for use with + :class:`~meshmode.discretization.Discretization`, + with one entry in the :class:`DOFArray` for each + :class:`~meshmode.discretization.ElementGroupBase`. + It is derived from :class:`numpy.ndarray` with dtype object ("an object array"). + The entries in this array are further arrays managed by :attr:`array_context`. + + The main purpose of this class is to better describe the data structure, + i.e. when a :class:`DOFArray` occurs inside of further numpy object array, + the level representing the array of element groups can be recognized (by + people and programs). + + .. attribute:: array_context + + An :class:`meshmode.array_context.ArrayContext`. + + .. attribute:: entry_dtype + + The (assumed uniform) :class:`numpy.dtype` of the group arrays + contained in this instance. + + .. automethod:: from_list + """ + + # Follows https://numpy.org/devdocs/user/basics.subclassing.html + + def __new__(cls, actx: Optional[ArrayContext], input_array): + if not (actx is None or isinstance(actx, ArrayContext)): + raise TypeError("actx must be of type ArrayContext") + + result = np.asarray(input_array).view(cls) + if len(result.shape) != 1: + raise ValueError("DOFArray instances must have one-dimensional " + "shape, with one entry per element group") + + result.array_context = actx + return result + + def __array_finalize__(self, obj): + if obj is None: + return + self.array_context = getattr(obj, "array_context", None) + + @property + def entry_dtype(self): + return single_valued(subary.dtype for subary in self.flat) + + @classmethod + def from_list(cls, actx: Optional[ArrayContext], res_list) -> "DOFArray": + r"""Create a :class:`DOFArray` from a list of arrays + (one per :class:`~meshmode.discretization.ElementGroupBase`). + + :arg actx: If *None*, the arrays in *res_list* must be + :meth:`~meshmode.array_context.ArrayContext.thaw`\ ed. + """ + if not (actx is None or isinstance(actx, ArrayContext)): + raise TypeError("actx must be of type ArrayContext") + + ngroups = len(res_list) + + result = np.empty(ngroups, dtype=object).view(cls) + result.array_context = actx + + # 'result[:] = res_list' may look tempting, however: + # https://github.com/numpy/numpy/issues/16564 + for igrp in range(ngroups): + result[igrp] = res_list[igrp] + + return result + +# }}} + + +def thaw(actx: ArrayContext, ary: np.ndarray) -> np.ndarray: + r"""Call :meth:`~meshmode.array_context.ArrayContext.thaw` on the element + group arrays making up the :class:`DOFArray`, using *actx*. + + Vectorizes over object arrays of :class:`DOFArray`\ s. + """ + if (isinstance(ary, np.ndarray) + and ary.dtype.char == "O" + and not isinstance(ary, DOFArray)): + return obj_array_vectorize(partial(thaw, actx), ary) + + if ary.array_context is not None: + raise ValueError("DOFArray passed to thaw is not frozen") + + return DOFArray.from_list(actx, [ + actx.thaw(subary) + for subary in ary + ]) + + +def freeze(ary: np.ndarray) -> np.ndarray: + r"""Call :meth:`~meshmode.array_context.ArrayContext.freeze` on the element + group arrays making up the :class:`DOFArray`, using the + :class:`~meshmode.array_context.ArrayContext` in *ary*. + + Vectorizes over object arrays of :class:`DOFArray`\ s. + """ + if (isinstance(ary, np.ndarray) + and ary.dtype.char == "O" + and not isinstance(ary, DOFArray)): + return obj_array_vectorize(freeze, ary) + + if ary.array_context is None: + raise ValueError("DOFArray passed to freeze is already frozen") + + return DOFArray.from_list(None, [ + ary.array_context.freeze(subary) + for subary in ary + ]) + + +def flatten(ary: np.ndarray) -> np.ndarray: + r"""Convert a :class:`DOFArray` into a "flat" array of degrees of freedom, + where the resulting type of the array is given by the + :attr:`DOFArray.array_context`. + + Array elements are laid out contiguously, with the element group + index varying slowest, element index next, and intra-element DOF + index fastest. + + Vectorizes over object arrays of :class:`DOFArray`\ s. + """ + if (isinstance(ary, np.ndarray) + and ary.dtype.char == "O" + and not isinstance(ary, DOFArray)): + return obj_array_vectorize(flatten, ary) + + group_sizes = [grp_ary.shape[0] * grp_ary.shape[1] for grp_ary in ary] + group_starts = np.cumsum([0] + group_sizes) + + actx = ary.array_context + + @memoize_in(actx, (flatten, "flatten_prg")) + def prg(): + return make_loopy_program( + "{[iel,idof]: 0<=iel np.ndarray: + r"""Convert a 'flat' array returned by :func:`flatten` back to a :class:`DOFArray`. + + Vectorizes over object arrays of :class:`DOFArray`\ s. + """ + if (isinstance(ary, np.ndarray) + and ary.dtype.char == "O" + and not isinstance(ary, DOFArray)): + return obj_array_vectorize( + lambda subary: unflatten(actx, discr, subary), + ary) + + @memoize_in(actx, (unflatten, "unflatten_prg")) + def prg(): + return make_loopy_program( + "{[iel,idof]: 0<=iel=2018.4", + "pytools>=2020.3", "pytest>=2.3", "loo.py>=2014.1", ], diff --git a/test/test_chained.py b/test/test_chained.py index 21ebbda4f202b5175f3d04e5c4eae2c95027b759..51c1e12c494f32af90e8a7fa98724188a04bc2ca 100644 --- a/test/test_chained.py +++ b/test/test_chained.py @@ -23,27 +23,25 @@ THE SOFTWARE. """ import numpy as np -import numpy.linalg as la import pyopencl as cl -import pyopencl.array # noqa -import pyopencl.clmath # noqa import pytest from pyopencl.tools import ( # noqa pytest_generate_tests_for_pyopencl as pytest_generate_tests) +from meshmode.array_context import PyOpenCLArrayContext +from meshmode.dof_array import thaw, flat_norm, flatten + import logging logger = logging.getLogger(__name__) -def create_discretization(queue, ndim, +def create_discretization(actx, ndim, nelements=42, mesh_name=None, order=4): - ctx = queue.context - # construct mesh if ndim == 2: from functools import partial @@ -80,13 +78,13 @@ def create_discretization(queue, ndim, from meshmode.discretization import Discretization from meshmode.discretization.poly_element import \ InterpolatoryQuadratureSimplexGroupFactory - discr = Discretization(ctx, mesh, + discr = Discretization(actx, mesh, InterpolatoryQuadratureSimplexGroupFactory(order)) return discr -def create_refined_connection(queue, discr, threshold=0.3): +def create_refined_connection(actx, discr, threshold=0.3): from meshmode.mesh.refinement import RefinerWithoutAdjacency from meshmode.discretization.connection import make_refinement_connection from meshmode.discretization.poly_element import \ @@ -97,20 +95,20 @@ def create_refined_connection(queue, discr, threshold=0.3): refiner.refine(flags) discr_order = discr.groups[0].order - connection = make_refinement_connection(refiner, discr, + connection = make_refinement_connection(actx, refiner, discr, InterpolatoryQuadratureSimplexGroupFactory(discr_order)) return connection -def create_face_connection(queue, discr): +def create_face_connection(actx, discr): from meshmode.discretization.connection import FACE_RESTR_ALL from meshmode.discretization.connection import make_face_restriction from meshmode.discretization.poly_element import \ InterpolatoryQuadratureSimplexGroupFactory discr_order = discr.groups[0].order - connection = make_face_restriction(discr, + connection = make_face_restriction(actx, discr, InterpolatoryQuadratureSimplexGroupFactory(discr_order), FACE_RESTR_ALL, per_face_groups=True) @@ -126,19 +124,20 @@ def test_chained_batch_table(ctx_factory, ndim, visualize=False): ctx = ctx_factory() queue = cl.CommandQueue(ctx) + actx = PyOpenCLArrayContext(queue) - discr = create_discretization(queue, ndim) + discr = create_discretization(actx, ndim) connections = [] - conn = create_refined_connection(queue, discr) + conn = create_refined_connection(actx, discr) connections.append(conn) - conn = create_refined_connection(queue, conn.to_discr) + conn = create_refined_connection(actx, conn.to_discr) connections.append(conn) from meshmode.discretization.connection import ChainedDiscretizationConnection chained = ChainedDiscretizationConnection(connections) conn = chained.connections[0] - el_table = _build_element_lookup_table(queue, conn) + el_table = _build_element_lookup_table(actx, conn) for igrp, grp in enumerate(conn.groups): for ibatch, batch in enumerate(grp.batches): ifrom = batch.from_element_indices.get(queue) @@ -156,15 +155,15 @@ def test_chained_new_group_table(ctx_factory, ndim, visualize=False): ctx = ctx_factory() queue = cl.CommandQueue(ctx) + actx = PyOpenCLArrayContext(queue) - discr = create_discretization(queue, ndim, + discr = create_discretization(actx, ndim, nelements=8, - mesh_order=2, - discr_order=2) + order=2) connections = [] - conn = create_refined_connection(queue, discr) + conn = create_refined_connection(actx, discr) connections.append(conn) - conn = create_refined_connection(queue, conn.to_discr) + conn = create_refined_connection(actx, conn.to_discr) connections.append(conn) grp_to_grp, grp_info = _build_new_group_table(connections[0], @@ -202,13 +201,13 @@ def test_chained_new_group_table(ctx_factory, ndim, visualize=False): def test_chained_connection(ctx_factory, ndim, visualize=False): ctx = ctx_factory() queue = cl.CommandQueue(ctx) + actx = PyOpenCLArrayContext(queue) - discr = create_discretization(queue, ndim, - nelements=10) + discr = create_discretization(actx, ndim, nelements=10) connections = [] - conn = create_refined_connection(queue, discr, threshold=np.inf) + conn = create_refined_connection(actx, discr, threshold=np.inf) connections.append(conn) - conn = create_refined_connection(queue, conn.to_discr, threshold=np.inf) + conn = create_refined_connection(actx, conn.to_discr, threshold=np.inf) connections.append(conn) from meshmode.discretization.connection import \ @@ -216,15 +215,15 @@ def test_chained_connection(ctx_factory, ndim, visualize=False): chained = ChainedDiscretizationConnection(connections) def f(x): - from six.moves import reduce - return 0.1 * reduce(lambda x, y: x * cl.clmath.sin(5 * y), x) + from functools import reduce + return 0.1 * reduce(lambda x, y: x * actx.np.sin(5 * y), x) - x = connections[0].from_discr.nodes().with_queue(queue) + x = thaw(actx, connections[0].from_discr.nodes()) fx = f(x) - f1 = chained(queue, fx).get(queue) - f2 = connections[1](queue, connections[0](queue, fx)).get(queue) + f1 = chained(fx) + f2 = connections[1](connections[0](fx)) - assert np.allclose(f1, f2) + assert flat_norm(f1-f2, np.inf) / flat_norm(f2) < 1e-11 @pytest.mark.skip(reason='slow test') @@ -235,12 +234,13 @@ def test_chained_full_resample_matrix(ctx_factory, ndim, visualize=False): ctx = ctx_factory() queue = cl.CommandQueue(ctx) + actx = PyOpenCLArrayContext(queue) - discr = create_discretization(queue, ndim) + discr = create_discretization(actx, ndim) connections = [] - conn = create_refined_connection(queue, discr) + conn = create_refined_connection(actx, discr) connections.append(conn) - conn = create_refined_connection(queue, conn.to_discr) + conn = create_refined_connection(actx, conn.to_discr) connections.append(conn) from meshmode.discretization.connection import \ @@ -249,15 +249,15 @@ def test_chained_full_resample_matrix(ctx_factory, ndim, visualize=False): def f(x): from six.moves import reduce - return 0.1 * reduce(lambda x, y: x * cl.clmath.sin(5 * y), x) + return 0.1 * reduce(lambda x, y: x * actx.np.sin(5 * y), x) - resample_mat = make_full_resample_matrix(queue, chained).get(queue) + resample_mat = actx.to_numpy(make_full_resample_matrix(actx, chained)) - x = connections[0].from_discr.nodes().with_queue(queue) + x = thaw(actx, connections[0].from_discr.nodes()) fx = f(x) - f1 = np.dot(resample_mat, fx.get(queue)) - f2 = chained(queue, fx).get(queue) - f3 = connections[1](queue, connections[0](queue, fx)).get(queue) + f1 = resample_mat @ actx.to_numpy(flatten(fx)) + f2 = actx.to_numpy(flatten(chained(fx))) + f3 = actx.to_numpy(flatten(connections[1](connections[0](fx)))) assert np.allclose(f1, f2) assert np.allclose(f2, f3) @@ -273,25 +273,26 @@ def test_chained_to_direct(ctx_factory, ndim, chain_type, ctx = ctx_factory() queue = cl.CommandQueue(ctx) + actx = PyOpenCLArrayContext(queue) - discr = create_discretization(queue, ndim, nelements=nelements) + discr = create_discretization(actx, ndim, nelements=nelements) connections = [] if chain_type == 1: - conn = create_refined_connection(queue, discr) + conn = create_refined_connection(actx, discr) connections.append(conn) - conn = create_refined_connection(queue, conn.to_discr) + conn = create_refined_connection(actx, conn.to_discr) connections.append(conn) elif chain_type == 2: - conn = create_refined_connection(queue, discr) + conn = create_refined_connection(actx, discr) connections.append(conn) - conn = create_refined_connection(queue, conn.to_discr) + conn = create_refined_connection(actx, conn.to_discr) connections.append(conn) - conn = create_refined_connection(queue, conn.to_discr) + conn = create_refined_connection(actx, conn.to_discr) connections.append(conn) elif chain_type == 3 and ndim == 3: - conn = create_refined_connection(queue, discr, threshold=np.inf) + conn = create_refined_connection(actx, discr, threshold=np.inf) connections.append(conn) - conn = create_face_connection(queue, conn.to_discr) + conn = create_face_connection(actx, conn.to_discr) connections.append(conn) else: raise ValueError('unknown test case') @@ -301,7 +302,7 @@ def test_chained_to_direct(ctx_factory, ndim, chain_type, chained = ChainedDiscretizationConnection(connections) t_start = time.time() - direct = flatten_chained_connection(queue, chained) + direct = flatten_chained_connection(actx, chained) t_end = time.time() if visualize: print('[TIME] Flatten: {:.5e}'.format(t_end - t_start)) @@ -317,19 +318,19 @@ def test_chained_to_direct(ctx_factory, ndim, chain_type, def f(x): from six.moves import reduce - return 0.1 * reduce(lambda x, y: x * cl.clmath.sin(5 * y), x) + return 0.1 * reduce(lambda x, y: x * actx.np.sin(5 * y), x) - x = connections[0].from_discr.nodes().with_queue(queue) + x = thaw(actx, connections[0].from_discr.nodes()) fx = f(x) t_start = time.time() - f1 = direct(queue, fx).get(queue) + f1 = actx.to_numpy(flatten(direct(fx))) t_end = time.time() if visualize: print('[TIME] Direct: {:.5e}'.format(t_end - t_start)) t_start = time.time() - f2 = chained(queue, fx).get(queue) + f2 = actx.to_numpy(flatten(chained(fx))) t_end = time.time() if visualize: print('[TIME] Chained: {:.5e}'.format(t_end - t_start)) @@ -351,27 +352,28 @@ def test_chained_to_direct(ctx_factory, ndim, chain_type, @pytest.mark.parametrize(("ndim", "mesh_name"), [ (2, "starfish"), (3, "torus")]) -def test_reversed_chained_connection(ctx_factory, ndim, mesh_name, visualize=False): +def test_reversed_chained_connection(ctx_factory, ndim, mesh_name): ctx = ctx_factory() queue = cl.CommandQueue(ctx) + actx = PyOpenCLArrayContext(queue) def run(nelements, order): - discr = create_discretization(queue, ndim, + discr = create_discretization(actx, ndim, nelements=nelements, order=order, mesh_name=mesh_name) threshold = 1.0 connections = [] - conn = create_refined_connection(queue, + conn = create_refined_connection(actx, discr, threshold=threshold) connections.append(conn) if ndim == 2: # NOTE: additional refinement makes the 3D meshes explode in size - conn = create_refined_connection(queue, + conn = create_refined_connection(actx, conn.to_discr, threshold=threshold) connections.append(conn) - conn = create_refined_connection(queue, + conn = create_refined_connection(actx, conn.to_discr, threshold=threshold) connections.append(conn) @@ -383,21 +385,19 @@ def test_reversed_chained_connection(ctx_factory, ndim, mesh_name, visualize=Fal reverse = L2ProjectionInverseDiscretizationConnection(chained) # create test vector - from_nodes = chained.from_discr.nodes().with_queue(queue) - to_nodes = chained.to_discr.nodes().with_queue(queue) + from_nodes = thaw(actx, chained.from_discr.nodes()) + to_nodes = thaw(actx, chained.to_discr.nodes()) from_x = 0 to_x = 0 for d in range(ndim): - from_x += cl.clmath.cos(from_nodes[d]) ** (d + 1) - to_x += cl.clmath.cos(to_nodes[d]) ** (d + 1) - - from_interp = reverse(queue, to_x) + from_x += actx.np.cos(from_nodes[d]) ** (d + 1) + to_x += actx.np.cos(to_nodes[d]) ** (d + 1) - from_interp = from_interp.get(queue) - from_x = from_x.get(queue) + from_interp = reverse(to_x) - return 1.0 / nelements, la.norm(from_interp - from_x) / la.norm(from_x) + return (1.0 / nelements, + flat_norm(from_interp - from_x, np.inf) / flat_norm(from_x, np.inf)) from pytools.convergence import EOCRecorder eoc = EOCRecorder() diff --git a/test/test_meshmode.py b/test/test_meshmode.py index 047362d20c1b297e22eed758fbe633bee01a82cc..a888beff918cc0ef1ccc002fa2f9878123999911 100644 --- a/test/test_meshmode.py +++ b/test/test_meshmode.py @@ -26,8 +26,6 @@ from six.moves import range import numpy as np import numpy.linalg as la import pyopencl as cl -import pyopencl.array # noqa -import pyopencl.clmath # noqa from pyopencl.tools import ( # noqa pytest_generate_tests_for_pyopencl @@ -40,6 +38,8 @@ from meshmode.discretization.poly_element import ( PolynomialEquidistantSimplexGroupFactory, ) from meshmode.mesh import Mesh, BTAG_ALL +from meshmode.array_context import PyOpenCLArrayContext +from meshmode.dof_array import thaw, flat_norm, flatten from meshmode.discretization.connection import \ FACE_RESTR_ALL, FACE_RESTR_INTERIOR import meshmode.mesh.generation as mgen @@ -59,7 +59,8 @@ def test_circle_mesh(do_plot=False): FileSource("circle.step"), 2, order=2, force_ambient_dim=2, other_options=[ - "-string", "Mesh.CharacteristicLengthMax = 0.05;"] + "-string", "Mesh.CharacteristicLengthMax = 0.05;"], + target_unit="MM", ) print("END GEN") print(mesh.nelements) @@ -223,6 +224,7 @@ def test_boundary_interpolation(ctx_factory, group_factory, boundary_tag, mesh_name, dim, mesh_pars, per_face_groups): cl_ctx = ctx_factory() queue = cl.CommandQueue(cl_ctx) + actx = PyOpenCLArrayContext(queue) from meshmode.discretization import Discretization from meshmode.discretization.connection import ( @@ -234,7 +236,7 @@ def test_boundary_interpolation(ctx_factory, group_factory, boundary_tag, order = 4 def f(x): - return 0.1*cl.clmath.sin(30*x) + return 0.1*actx.np.sin(30*x) for mesh_par in mesh_pars: # {{{ get mesh @@ -267,32 +269,31 @@ def test_boundary_interpolation(ctx_factory, group_factory, boundary_tag, # }}} - vol_discr = Discretization(cl_ctx, mesh, - group_factory(order)) + vol_discr = Discretization(actx, mesh, group_factory(order)) print("h=%s -> %d elements" % ( h, sum(mgrp.nelements for mgrp in mesh.groups))) - x = vol_discr.nodes()[0].with_queue(queue) + x = thaw(actx, vol_discr.nodes()[0]) vol_f = f(x) bdry_connection = make_face_restriction( - vol_discr, group_factory(order), + actx, vol_discr, group_factory(order), boundary_tag, per_face_groups=per_face_groups) - check_connection(bdry_connection) + check_connection(actx, bdry_connection) bdry_discr = bdry_connection.to_discr - bdry_x = bdry_discr.nodes()[0].with_queue(queue) + bdry_x = thaw(actx, bdry_discr.nodes()[0]) bdry_f = f(bdry_x) - bdry_f_2 = bdry_connection(queue, vol_f) + bdry_f_2 = bdry_connection(vol_f) if mesh_name == "blob" and dim == 2 and mesh.nelements < 500: - mat = bdry_connection.full_resample_matrix(queue).get(queue) - bdry_f_2_by_mat = mat.dot(vol_f.get()) + mat = actx.to_numpy(bdry_connection.full_resample_matrix(actx)) + bdry_f_2_by_mat = mat.dot(actx.to_numpy(flatten(vol_f))) - mat_error = la.norm(bdry_f_2.get(queue=queue) - bdry_f_2_by_mat) + mat_error = la.norm(actx.to_numpy(flatten(bdry_f_2)) - bdry_f_2_by_mat) assert mat_error < 1e-14, mat_error - err = la.norm((bdry_f-bdry_f_2).get(), np.inf) + err = flat_norm(bdry_f-bdry_f_2, np.inf) eoc_rec.add_data_point(h, err) order_slack = 0.75 if mesh_name == "blob" else 0.5 @@ -316,6 +317,7 @@ def test_all_faces_interpolation(ctx_factory, mesh_name, dim, mesh_pars, per_face_groups): cl_ctx = ctx_factory() queue = cl.CommandQueue(cl_ctx) + actx = PyOpenCLArrayContext(queue) from meshmode.discretization import Discretization from meshmode.discretization.connection import ( @@ -328,7 +330,7 @@ def test_all_faces_interpolation(ctx_factory, mesh_name, dim, mesh_pars, order = 4 def f(x): - return 0.1*cl.clmath.sin(30*x) + return 0.1*actx.np.sin(30*x) for mesh_par in mesh_pars: # {{{ get mesh @@ -344,7 +346,8 @@ def test_all_faces_interpolation(ctx_factory, mesh_name, dim, mesh_pars, FileSource("blob-2d.step"), 2, order=order, force_ambient_dim=2, other_options=[ - "-string", "Mesh.CharacteristicLengthMax = %s;" % h] + "-string", "Mesh.CharacteristicLengthMax = %s;" % h], + target_unit="MM", ) print("END GEN") elif mesh_name == "warp": @@ -357,20 +360,20 @@ def test_all_faces_interpolation(ctx_factory, mesh_name, dim, mesh_pars, # }}} - vol_discr = Discretization(cl_ctx, mesh, + vol_discr = Discretization(actx, mesh, PolynomialWarpAndBlendGroupFactory(order)) print("h=%s -> %d elements" % ( h, sum(mgrp.nelements for mgrp in mesh.groups))) all_face_bdry_connection = make_face_restriction( - vol_discr, PolynomialWarpAndBlendGroupFactory(order), + actx, vol_discr, PolynomialWarpAndBlendGroupFactory(order), FACE_RESTR_ALL, per_face_groups=per_face_groups) all_face_bdry_discr = all_face_bdry_connection.to_discr for ito_grp, ceg in enumerate(all_face_bdry_connection.groups): for ibatch, batch in enumerate(ceg.batches): assert np.array_equal( - batch.from_element_indices.get(queue), + actx.to_numpy(actx.thaw(batch.from_element_indices)), np.arange(vol_discr.mesh.nelements)) if per_face_groups: @@ -378,31 +381,31 @@ def test_all_faces_interpolation(ctx_factory, mesh_name, dim, mesh_pars, else: assert ibatch == batch.to_element_face - all_face_x = all_face_bdry_discr.nodes()[0].with_queue(queue) + all_face_x = thaw(actx, all_face_bdry_discr.nodes()[0]) all_face_f = f(all_face_x) - all_face_f_2 = all_face_bdry_discr.zeros(queue) + all_face_f_2 = all_face_bdry_discr.zeros(actx) for boundary_tag in [ BTAG_ALL, FACE_RESTR_INTERIOR, ]: bdry_connection = make_face_restriction( - vol_discr, PolynomialWarpAndBlendGroupFactory(order), + actx, vol_discr, PolynomialWarpAndBlendGroupFactory(order), boundary_tag, per_face_groups=per_face_groups) bdry_discr = bdry_connection.to_discr - bdry_x = bdry_discr.nodes()[0].with_queue(queue) + bdry_x = thaw(actx, bdry_discr.nodes()[0]) bdry_f = f(bdry_x) all_face_embedding = make_face_to_all_faces_embedding( - bdry_connection, all_face_bdry_discr) + actx, bdry_connection, all_face_bdry_discr) - check_connection(all_face_embedding) + check_connection(actx, all_face_embedding) - all_face_f_2 += all_face_embedding(queue, bdry_f) + all_face_f_2 = all_face_f_2 + all_face_embedding(bdry_f) - err = la.norm((all_face_f-all_face_f_2).get(), np.inf) + err = flat_norm(all_face_f-all_face_f_2, np.inf) eoc_rec.add_data_point(h, err) print(eoc_rec) @@ -431,6 +434,7 @@ def test_opposite_face_interpolation(ctx_factory, group_factory, cl_ctx = ctx_factory() queue = cl.CommandQueue(cl_ctx) + actx = PyOpenCLArrayContext(queue) from meshmode.discretization import Discretization from meshmode.discretization.connection import ( @@ -443,7 +447,7 @@ def test_opposite_face_interpolation(ctx_factory, group_factory, order = 5 def f(x): - return 0.1*cl.clmath.sin(30*x) + return 0.1*actx.np.sin(30*x) for mesh_par in mesh_pars: # {{{ get mesh @@ -467,7 +471,8 @@ def test_opposite_face_interpolation(ctx_factory, group_factory, FileSource("blob-2d.step"), 2, order=order, force_ambient_dim=2, other_options=[ - "-string", "Mesh.CharacteristicLengthMax = %s;" % h] + "-string", "Mesh.CharacteristicLengthMax = %s;" % h], + target_unit="MM", ) print("END GEN") elif mesh_name == "warp": @@ -480,24 +485,24 @@ def test_opposite_face_interpolation(ctx_factory, group_factory, # }}} - vol_discr = Discretization(cl_ctx, mesh, + vol_discr = Discretization(actx, mesh, group_factory(order)) print("h=%s -> %d elements" % ( h, sum(mgrp.nelements for mgrp in mesh.groups))) bdry_connection = make_face_restriction( - vol_discr, group_factory(order), + actx, vol_discr, group_factory(order), FACE_RESTR_INTERIOR) bdry_discr = bdry_connection.to_discr - opp_face = make_opposite_face_connection(bdry_connection) - check_connection(opp_face) + opp_face = make_opposite_face_connection(actx, bdry_connection) + check_connection(actx, opp_face) - bdry_x = bdry_discr.nodes()[0].with_queue(queue) + bdry_x = thaw(actx, bdry_discr.nodes()[0]) bdry_f = f(bdry_x) - bdry_f_2 = opp_face(queue, bdry_f) + bdry_f_2 = opp_face(bdry_f) - err = la.norm((bdry_f-bdry_f_2).get(), np.inf) + err = flat_norm(bdry_f-bdry_f_2, np.inf) eoc_rec.add_data_point(h, err) print(eoc_rec) @@ -518,7 +523,8 @@ def test_element_orientation(): mesh = generate_gmsh( FileSource("blob-2d.step"), 2, order=mesh_order, force_ambient_dim=2, - other_options=["-string", "Mesh.CharacteristicLengthMax = 0.02;"] + other_options=["-string", "Mesh.CharacteristicLengthMax = 0.02;"], + target_unit="MM", ) from meshmode.mesh.processing import (perform_flips, @@ -555,13 +561,14 @@ def test_3d_orientation(ctx_factory, what, mesh_gen_func, visualize=False): ctx = ctx_factory() queue = cl.CommandQueue(ctx) + actx = PyOpenCLArrayContext(queue) mesh = mesh_gen_func() logger.info("%d elements" % mesh.nelements) from meshmode.discretization import Discretization - discr = Discretization(ctx, mesh, + discr = Discretization(actx, mesh, PolynomialWarpAndBlendGroupFactory(1)) from pytential import bind, sym @@ -582,17 +589,18 @@ def test_3d_orientation(ctx_factory, what, mesh_gen_func, visualize=False): normal_outward_expr = ( sym.normal(mesh.ambient_dim) | sym.nodes(mesh.ambient_dim)) - normal_outward_check = bind(discr, normal_outward_expr)(queue).as_scalar() > 0 + normal_outward_check = actx.to_numpy( + flatten(bind(discr, normal_outward_expr)(actx).as_scalar())) > 0 - assert normal_outward_check.get().all(), normal_outward_check.get() + assert normal_outward_check.all(), normal_outward_check # }}} - normals = bind(discr, sym.normal(mesh.ambient_dim).xproject(1))(queue) + normals = bind(discr, sym.normal(mesh.ambient_dim).xproject(1))(actx) if visualize: from meshmode.discretization.visualization import make_visualizer - vis = make_visualizer(queue, discr, 1) + vis = make_visualizer(actx, discr, 1) vis.write_vtk_file("normals.vtu", [ ("normals", normals), @@ -616,7 +624,8 @@ def test_merge_and_map(ctx_factory, visualize=False): mesh = generate_gmsh( FileSource("blob-2d.step"), 2, order=mesh_order, force_ambient_dim=2, - other_options=["-string", "Mesh.CharacteristicLengthMax = 0.02;"] + other_options=["-string", "Mesh.CharacteristicLengthMax = 0.02;"], + target_unit="MM", ) discr_grp_factory = PolynomialWarpAndBlendGroupFactory(3) @@ -664,6 +673,7 @@ def test_sanity_single_element(ctx_factory, dim, order, visualize=False): cl_ctx = ctx_factory() queue = cl.CommandQueue(cl_ctx) + actx = PyOpenCLArrayContext(queue) from modepy.tools import unit_vertices vertices = unit_vertices(dim).T.copy() @@ -684,21 +694,20 @@ def test_sanity_single_element(ctx_factory, dim, order, visualize=False): from meshmode.discretization import Discretization from meshmode.discretization.poly_element import \ PolynomialWarpAndBlendGroupFactory - vol_discr = Discretization(cl_ctx, mesh, + vol_discr = Discretization(actx, mesh, PolynomialWarpAndBlendGroupFactory(order+3)) # {{{ volume calculation check - vol_x = vol_discr.nodes().with_queue(queue) + vol_x = thaw(actx, vol_discr.nodes()) - vol_one = vol_x[0].copy() - vol_one.fill(1) + vol_one = vol_x[0] * 0 + 1 from pytential import norm, integral # noqa from pytools import factorial true_vol = 1/factorial(dim) * 2**dim - comp_vol = integral(vol_discr, queue, vol_one) + comp_vol = integral(vol_discr, vol_one) rel_vol_err = abs(true_vol - comp_vol) / true_vol assert rel_vol_err < 1e-12 @@ -709,7 +718,7 @@ def test_sanity_single_element(ctx_factory, dim, order, visualize=False): from meshmode.discretization.connection import make_face_restriction bdry_connection = make_face_restriction( - vol_discr, PolynomialWarpAndBlendGroupFactory(order + 3), + actx, vol_discr, PolynomialWarpAndBlendGroupFactory(order + 3), BTAG_ALL) bdry_discr = bdry_connection.to_discr @@ -719,12 +728,12 @@ def test_sanity_single_element(ctx_factory, dim, order, visualize=False): from meshmode.discretization.visualization import make_visualizer #vol_vis = make_visualizer(queue, vol_discr, 4) - bdry_vis = make_visualizer(queue, bdry_discr, 4) + bdry_vis = make_visualizer(actx, bdry_discr, 4) # }}} from pytential import bind, sym - bdry_normals = bind(bdry_discr, sym.normal(dim))(queue).as_vector(dtype=object) + bdry_normals = bind(bdry_discr, sym.normal(dim))(actx).as_vector(dtype=object) if visualize: bdry_vis.write_vtk_file("boundary.vtu", [ @@ -734,9 +743,10 @@ def test_sanity_single_element(ctx_factory, dim, order, visualize=False): normal_outward_check = bind(bdry_discr, sym.normal(dim) | (sym.nodes(dim) + 0.5*sym.ones_vec(dim)), - )(queue).as_scalar() > 0 + )(actx).as_scalar() - assert normal_outward_check.get().all(), normal_outward_check.get() + normal_outward_check = actx.to_numpy(flatten(normal_outward_check) > 0) + assert normal_outward_check.all(), normal_outward_check # }}} @@ -752,6 +762,7 @@ def test_sanity_qhull_nd(ctx_factory, dim, order): ctx = ctx_factory() queue = cl.CommandQueue(ctx) + actx = PyOpenCLArrayContext(queue) from scipy.spatial import Delaunay verts = np.random.rand(1000, dim) @@ -762,28 +773,28 @@ def test_sanity_qhull_nd(ctx_factory, dim, order): fix_orientation=True) from meshmode.discretization import Discretization - low_discr = Discretization(ctx, mesh, + low_discr = Discretization(actx, mesh, PolynomialEquidistantSimplexGroupFactory(order)) - high_discr = Discretization(ctx, mesh, + high_discr = Discretization(actx, mesh, PolynomialEquidistantSimplexGroupFactory(order+1)) from meshmode.discretization.connection import make_same_mesh_connection - cnx = make_same_mesh_connection(high_discr, low_discr) + cnx = make_same_mesh_connection(actx, high_discr, low_discr) def f(x): - return 0.1*cl.clmath.sin(x) + return 0.1*actx.np.sin(x) - x_low = low_discr.nodes()[0].with_queue(queue) + x_low = thaw(actx, low_discr.nodes()[0]) f_low = f(x_low) - x_high = high_discr.nodes()[0].with_queue(queue) + x_high = thaw(actx, high_discr.nodes()[0]) f_high_ref = f(x_high) - f_high_num = cnx(queue, f_low) + f_high_num = cnx(f_low) - err = (f_high_ref-f_high_num).get() - - err = la.norm(err, np.inf)/la.norm(f_high_ref.get(), np.inf) + err = ( + flat_norm(f_high_ref-f_high_num, np.inf) + / flat_norm(f_high_ref, np.inf)) print(err) assert err < 1e-2 @@ -799,14 +810,14 @@ def test_sanity_qhull_nd(ctx_factory, dim, order): ("ball-radius-1.step", 3), ]) @pytest.mark.parametrize("mesh_order", [1, 2]) -def test_sanity_balls(ctx_factory, src_file, dim, mesh_order, - visualize=False): +def test_sanity_balls(ctx_factory, src_file, dim, mesh_order, visualize=False): pytest.importorskip("pytential") logging.basicConfig(level=logging.INFO) ctx = ctx_factory() queue = cl.CommandQueue(ctx) + actx = PyOpenCLArrayContext(queue) from pytools.convergence import EOCRecorder vol_eoc_rec = EOCRecorder() @@ -822,18 +833,20 @@ def test_sanity_balls(ctx_factory, src_file, dim, mesh_order, mesh = generate_gmsh( FileSource(src_file), dim, order=mesh_order, other_options=["-string", "Mesh.CharacteristicLengthMax = %g;" % h], - force_ambient_dim=dim) + force_ambient_dim=dim, + target_unit="MM") logger.info("%d elements" % mesh.nelements) # {{{ discretizations and connections from meshmode.discretization import Discretization - vol_discr = Discretization(ctx, mesh, + vol_discr = Discretization(actx, mesh, InterpolatoryQuadratureSimplexGroupFactory(quad_order)) from meshmode.discretization.connection import make_face_restriction bdry_connection = make_face_restriction( + actx, vol_discr, InterpolatoryQuadratureSimplexGroupFactory(quad_order), BTAG_ALL) @@ -844,8 +857,8 @@ def test_sanity_balls(ctx_factory, src_file, dim, mesh_order, # {{{ visualizers from meshmode.discretization.visualization import make_visualizer - vol_vis = make_visualizer(queue, vol_discr, 20) - bdry_vis = make_visualizer(queue, bdry_discr, 20) + vol_vis = make_visualizer(actx, vol_discr, 20) + bdry_vis = make_visualizer(actx, bdry_discr, 20) # }}} @@ -853,27 +866,25 @@ def test_sanity_balls(ctx_factory, src_file, dim, mesh_order, true_surf = 2*np.pi**(dim/2)/gamma(dim/2) true_vol = true_surf/dim - vol_x = vol_discr.nodes().with_queue(queue) + vol_x = thaw(actx, vol_discr.nodes()) - vol_one = vol_x[0].copy() - vol_one.fill(1) + vol_one = vol_x[0]*0 + 1 from pytential import norm, integral # noqa - comp_vol = integral(vol_discr, queue, vol_one) + comp_vol = integral(vol_discr, vol_one) rel_vol_err = abs(true_vol - comp_vol) / true_vol vol_eoc_rec.add_data_point(h, rel_vol_err) print("VOL", true_vol, comp_vol) - bdry_x = bdry_discr.nodes().with_queue(queue) + bdry_x = thaw(actx, bdry_discr.nodes()) - bdry_one_exact = bdry_x[0].copy() - bdry_one_exact.fill(1) + bdry_one_exact = bdry_x[0] * 0 + 1 - bdry_one = bdry_connection(queue, vol_one).with_queue(queue) - intp_err = norm(bdry_discr, queue, bdry_one-bdry_one_exact) + bdry_one = bdry_connection(vol_one) + intp_err = norm(bdry_discr, bdry_one-bdry_one_exact) assert intp_err < 1e-14 - comp_surf = integral(bdry_discr, queue, bdry_one) + comp_surf = integral(bdry_discr, bdry_one) rel_surf_err = abs(true_surf - comp_surf) / true_surf surf_eoc_rec.add_data_point(h, rel_surf_err) print("SURF", true_surf, comp_surf) @@ -884,7 +895,7 @@ def test_sanity_balls(ctx_factory, src_file, dim, mesh_order, ("area_el", bind( vol_discr, sym.area_element(mesh.ambient_dim, mesh.ambient_dim)) - (queue)), + (actx)), ]) bdry_vis.write_vtk_file("boundary-h=%g.vtu" % h, [("f", bdry_one)]) @@ -892,9 +903,10 @@ def test_sanity_balls(ctx_factory, src_file, dim, mesh_order, normal_outward_check = bind(bdry_discr, sym.normal(mesh.ambient_dim) | sym.nodes(mesh.ambient_dim), - )(queue).as_scalar() > 0 + )(actx).as_scalar() - assert normal_outward_check.get().all(), normal_outward_check.get() + normal_outward_check = actx.to_numpy(flatten(normal_outward_check) > 0) + assert normal_outward_check.all(), normal_outward_check # }}} @@ -1044,6 +1056,7 @@ def test_quad_mesh_2d(): """, ["blob-2d.step"]), force_ambient_dim=2, + target_unit="MM", ) print("END GEN") print(mesh.nelements) @@ -1135,7 +1148,7 @@ def test_quad_multi_element(): # {{{ test_vtk_overwrite -def test_vtk_overwrite(ctx_getter): +def test_vtk_overwrite(ctx_factory): pytest.importorskip("pyvisfile") def _try_write_vtk(writer, obj): @@ -1154,8 +1167,9 @@ def test_vtk_overwrite(ctx_getter): if os.path.exists(filename): os.remove(filename) - ctx = ctx_getter() + ctx = ctx_factory() queue = cl.CommandQueue(ctx) + actx = PyOpenCLArrayContext(queue) target_order = 7 @@ -1166,7 +1180,7 @@ def test_vtk_overwrite(ctx_getter): from meshmode.discretization.poly_element import \ InterpolatoryQuadratureSimplexGroupFactory discr = Discretization( - queue.context, mesh, + actx, mesh, InterpolatoryQuadratureSimplexGroupFactory(target_order)) from meshmode.discretization.visualization import make_visualizer @@ -1174,7 +1188,7 @@ def test_vtk_overwrite(ctx_getter): write_nodal_adjacency_vtk_file from meshmode.mesh.visualization import write_vertex_vtk_file - vis = make_visualizer(queue, discr, 1) + vis = make_visualizer(actx, discr, 1) _try_write_vtk(vis.write_vtk_file, discr) _try_write_vtk(lambda x, y, **kwargs: @@ -1197,7 +1211,8 @@ def test_mesh_to_tikz(): FileSource("../test/blob-2d.step"), 2, order=order, force_ambient_dim=2, other_options=[ - "-string", "Mesh.CharacteristicLengthMax = %s;" % h] + "-string", "Mesh.CharacteristicLengthMax = %s;" % h], + target_unit="MM", ) from meshmode.mesh.visualization import mesh_to_tikz @@ -1227,6 +1242,7 @@ def test_affine_map(): def test_mesh_without_vertices(ctx_factory): ctx = ctx_factory() queue = cl.CommandQueue(ctx) + actx = PyOpenCLArrayContext(queue) # create a mesh from meshmode.mesh.generation import generate_icosphere @@ -1246,11 +1262,11 @@ def test_mesh_without_vertices(ctx_factory): from meshmode.discretization import Discretization from meshmode.discretization.poly_element import \ InterpolatoryQuadratureSimplexGroupFactory as GroupFactory - discr = Discretization(ctx, mesh, GroupFactory(4)) - discr.nodes().with_queue(queue) + discr = Discretization(actx, mesh, GroupFactory(4)) + thaw(actx, discr.nodes()) from meshmode.discretization.visualization import make_visualizer - make_visualizer(queue, discr, 4) + make_visualizer(actx, discr, 4) @pytest.mark.parametrize("curve_name", ["ellipse", "arc"]) @@ -1352,6 +1368,7 @@ def test_is_affine_group_check(mesh_name): def test_mesh_multiple_groups(ctx_factory, ambient_dim, visualize=False): ctx = ctx_factory() queue = cl.CommandQueue(ctx) + actx = PyOpenCLArrayContext(queue) order = 4 @@ -1385,16 +1402,16 @@ def test_mesh_multiple_groups(ctx_factory, ambient_dim, visualize=False): from meshmode.discretization import Discretization from meshmode.discretization.poly_element import \ PolynomialWarpAndBlendGroupFactory as GroupFactory - discr = Discretization(ctx, mesh, GroupFactory(order)) + discr = Discretization(actx, mesh, GroupFactory(order)) if visualize: - group_id = discr.empty(queue, dtype=np.int) + group_id = discr.empty(actx, dtype=np.int) for igrp, grp in enumerate(discr.groups): group_id_view = grp.view(group_id) group_id_view.fill(igrp) from meshmode.discretization.visualization import make_visualizer - vis = make_visualizer(queue, discr, vis_order=order) + vis = make_visualizer(actx, discr, vis_order=order) vis.write_vtk_file("test_mesh_multiple_groups.vtu", [ ("group_id", group_id) ], overwrite=True) @@ -1406,28 +1423,27 @@ def test_mesh_multiple_groups(ctx_factory, ambient_dim, visualize=False): make_opposite_face_connection, check_connection) for boundary_tag in [BTAG_ALL, FACE_RESTR_INTERIOR, FACE_RESTR_ALL]: - conn = make_face_restriction(discr, GroupFactory(order), + conn = make_face_restriction(actx, discr, GroupFactory(order), boundary_tag=boundary_tag, per_face_groups=False) - check_connection(conn) + check_connection(actx, conn) - bdry_f = conn.to_discr.empty(queue) - bdry_f.fill(1.0) + bdry_f = conn.to_discr.zeros(actx) + 1 if boundary_tag == FACE_RESTR_INTERIOR: - opposite = make_opposite_face_connection(conn) - check_connection(opposite) + opposite = make_opposite_face_connection(actx, conn) + check_connection(actx, opposite) - op_bdry_f = opposite(queue, bdry_f) - error = abs(cl.array.sum(bdry_f - op_bdry_f).get(queue)) + op_bdry_f = opposite(bdry_f) + error = flat_norm(bdry_f - op_bdry_f, np.inf) assert error < 1.0e-11, error if boundary_tag == FACE_RESTR_ALL: - embedding = make_face_to_all_faces_embedding(conn, conn.to_discr) - check_connection(embedding) + embedding = make_face_to_all_faces_embedding(actx, conn, conn.to_discr) + check_connection(actx, embedding) - em_bdry_f = embedding(queue, bdry_f) - error = abs(cl.array.sum(bdry_f - em_bdry_f).get(queue)) + em_bdry_f = embedding(bdry_f) + error = flat_norm(bdry_f - em_bdry_f) assert error < 1.0e-11, error diff --git a/test/test_partition.py b/test/test_partition.py index 54b9d88923249b553c334a278744e2f51ab14c44..1e8e6c7fb5d556807375dc3e132b4e1b1c3f3a09 100644 --- a/test/test_partition.py +++ b/test/test_partition.py @@ -27,10 +27,10 @@ THE SOFTWARE. from six.moves import range import numpy as np -import numpy.linalg as la import pyopencl as cl -import pyopencl.array # noqa -import pyopencl.clmath # noqa + +from meshmode.array_context import PyOpenCLArrayContext +from meshmode.dof_array import thaw, flatten, unflatten, flat_norm from pyopencl.tools import ( # noqa pytest_generate_tests_for_pyopencl @@ -69,6 +69,8 @@ def test_partition_interpolation(ctx_factory, dim, mesh_pars, group_factory = PolynomialWarpAndBlendGroupFactory cl_ctx = ctx_factory() queue = cl.CommandQueue(cl_ctx) + actx = PyOpenCLArrayContext(queue) + order = 4 from pytools.convergence import EOCRecorder @@ -80,7 +82,7 @@ def test_partition_interpolation(ctx_factory, dim, mesh_pars, eoc_rec[i, j] = EOCRecorder() def f(x): - return 10.*cl.clmath.sin(50.*x) + return 10.*actx.np.sin(50.*x) for n in mesh_pars: from meshmode.mesh.generation import generate_warped_rect_mesh @@ -107,7 +109,7 @@ def test_partition_interpolation(ctx_factory, dim, mesh_pars, partition_mesh(mesh, part_per_element, i)[0] for i in range(num_parts)] from meshmode.discretization import Discretization - vol_discrs = [Discretization(cl_ctx, part_meshes[i], group_factory(order)) + vol_discrs = [Discretization(actx, part_meshes[i], group_factory(order)) for i in range(num_parts)] from meshmode.mesh import BTAG_PARTITION @@ -120,23 +122,26 @@ def test_partition_interpolation(ctx_factory, dim, mesh_pars, continue # Mark faces within local_mesh that are connected to remote_mesh - local_bdry_conn = make_face_restriction(vol_discrs[i_local_part], + local_bdry_conn = make_face_restriction(actx, vol_discrs[i_local_part], group_factory(order), BTAG_PARTITION(i_remote_part)) # If these parts are not connected, don't bother checking the error - bdry_nodes = local_bdry_conn.to_discr.nodes() - if bdry_nodes.size == 0: + bdry_nelements = sum( + grp.nelements for grp in local_bdry_conn.to_discr.groups) + if bdry_nelements == 0: eoc_rec[i_local_part, i_remote_part] = None continue # Mark faces within remote_mesh that are connected to local_mesh - remote_bdry_conn = make_face_restriction(vol_discrs[i_remote_part], + remote_bdry_conn = make_face_restriction(actx, vol_discrs[i_remote_part], group_factory(order), BTAG_PARTITION(i_local_part)) - assert bdry_nodes.size == remote_bdry_conn.to_discr.nodes().size, \ - "partitions do not have the same number of connected nodes" + remote_bdry_nelements = sum( + grp.nelements for grp in remote_bdry_conn.to_discr.groups) + assert bdry_nelements == remote_bdry_nelements, \ + "partitions do not have the same number of connected elements" # Gather just enough information for the connection local_bdry = local_bdry_conn.to_discr @@ -166,27 +171,25 @@ def test_partition_interpolation(ctx_factory, dim, mesh_pars, for grp_batches in remote_batches] # Connect from remote_mesh to local_mesh - remote_to_local_conn = make_partition_connection(local_bdry_conn, - i_local_part, - remote_bdry, - remote_adj_groups, - remote_from_elem_faces, - remote_from_elem_indices) + remote_to_local_conn = make_partition_connection( + actx, local_bdry_conn, i_local_part, remote_bdry, + remote_adj_groups, remote_from_elem_faces, + remote_from_elem_indices) + # Connect from local mesh to remote mesh - local_to_remote_conn = make_partition_connection(remote_bdry_conn, - i_remote_part, - local_bdry, - local_adj_groups, - local_from_elem_faces, - local_from_elem_indices) - check_connection(remote_to_local_conn) - check_connection(local_to_remote_conn) - - true_local_points = f(local_bdry.nodes()[0].with_queue(queue)) - remote_points = local_to_remote_conn(queue, true_local_points) - local_points = remote_to_local_conn(queue, remote_points) - - err = la.norm((true_local_points - local_points).get(), np.inf) + local_to_remote_conn = make_partition_connection( + actx, remote_bdry_conn, i_remote_part, local_bdry, + local_adj_groups, local_from_elem_faces, + local_from_elem_indices) + + check_connection(actx, remote_to_local_conn) + check_connection(actx, local_to_remote_conn) + + true_local_points = f(thaw(actx, local_bdry.nodes()[0])) + remote_points = local_to_remote_conn(true_local_points) + local_points = remote_to_local_conn(remote_points) + + err = flat_norm(true_local_points - local_points, np.inf) eoc_rec[i_local_part, i_remote_part].add_data_point(1./n, err) for (i, j), e in eoc_rec.items(): @@ -359,9 +362,10 @@ def _test_mpi_boundary_swap(dim, order, num_groups): group_factory = PolynomialWarpAndBlendGroupFactory(order) cl_ctx = cl.create_some_context() queue = cl.CommandQueue(cl_ctx) + actx = PyOpenCLArrayContext(queue) from meshmode.discretization import Discretization - vol_discr = Discretization(cl_ctx, local_mesh, group_factory) + vol_discr = Discretization(actx, local_mesh, group_factory) from meshmode.distributed import get_connected_partitions connected_parts = get_connected_partitions(local_mesh) @@ -373,11 +377,11 @@ def _test_mpi_boundary_swap(dim, order, num_groups): from meshmode.mesh import BTAG_PARTITION for i_remote_part in connected_parts: local_bdry_conns[i_remote_part] = make_face_restriction( - vol_discr, group_factory, BTAG_PARTITION(i_remote_part)) + actx, vol_discr, group_factory, BTAG_PARTITION(i_remote_part)) setup_helper = bdry_setup_helpers[i_remote_part] = \ MPIBoundaryCommSetupHelper( - mpi_comm, queue, local_bdry_conns[i_remote_part], + mpi_comm, actx, local_bdry_conns[i_remote_part], i_remote_part, bdry_grp_factory=group_factory) setup_helper.post_sends() @@ -389,14 +393,14 @@ def _test_mpi_boundary_swap(dim, order, num_groups): if setup_helper.is_setup_ready(): assert bdry_setup_helpers.pop(i_remote_part) is setup_helper conn = setup_helper.complete_setup() - check_connection(conn) + check_connection(actx, conn) remote_to_local_bdry_conns[i_remote_part] = conn break # FIXME: Not ideal, busy-waits _test_data_transfer(mpi_comm, - queue, + actx, local_bdry_conns, remote_to_local_bdry_conns, connected_parts) @@ -405,12 +409,12 @@ def _test_mpi_boundary_swap(dim, order, num_groups): # TODO -def _test_data_transfer(mpi_comm, queue, local_bdry_conns, +def _test_data_transfer(mpi_comm, actx, local_bdry_conns, remote_to_local_bdry_conns, connected_parts): from mpi4py import MPI def f(x): - return 10*cl.clmath.sin(20.*x) + return 10*actx.np.sin(20.*x) ''' Here is a simplified example of what happens from @@ -435,13 +439,13 @@ def _test_data_transfer(mpi_comm, queue, local_bdry_conns, for i_remote_part in connected_parts: conn = remote_to_local_bdry_conns[i_remote_part] bdry_discr = local_bdry_conns[i_remote_part].to_discr - bdry_x = bdry_discr.nodes()[0].with_queue(queue=queue) + bdry_x = thaw(actx, bdry_discr.nodes()[0]) true_local_f = f(bdry_x) - remote_f = conn(queue, true_local_f) + remote_f = conn(true_local_f) # 2. - send_reqs.append(mpi_comm.isend(remote_f.get(queue=queue), + send_reqs.append(mpi_comm.isend(actx.to_numpy(flatten(remote_f)), dest=i_remote_part, tag=TAG_SEND_REMOTE_NODES)) @@ -471,12 +475,9 @@ def _test_data_transfer(mpi_comm, queue, local_bdry_conns, send_reqs = [] for i_remote_part in connected_parts: conn = remote_to_local_bdry_conns[i_remote_part] - local_f_np = remote_to_local_f_data[i_remote_part] - local_f_cl = cl.array.Array(queue, - shape=local_f_np.shape, - dtype=local_f_np.dtype) - local_f_cl.set(local_f_np) - remote_f = conn(queue, local_f_cl).get(queue=queue) + local_f = unflatten(actx, conn.from_discr, + actx.from_numpy(remote_to_local_f_data[i_remote_part])) + remote_f = actx.to_numpy(flatten(conn(local_f))) # 5. send_reqs.append(mpi_comm.isend(remote_f, @@ -508,9 +509,9 @@ def _test_data_transfer(mpi_comm, queue, local_bdry_conns, # 7. for i_remote_part in connected_parts: bdry_discr = local_bdry_conns[i_remote_part].to_discr - bdry_x = bdry_discr.nodes()[0].with_queue(queue=queue) + bdry_x = thaw(actx, bdry_discr.nodes()[0]) - true_local_f = f(bdry_x).get(queue=queue) + true_local_f = actx.to_numpy(flatten(f(bdry_x))) local_f = local_f_data[i_remote_part] from numpy.linalg import norm @@ -554,7 +555,7 @@ if __name__ == "__main__": if len(sys.argv) > 1: exec(sys.argv[1]) else: - from py.test.cmdline import main + from pytest import main main([__file__]) # vim: fdm=marker diff --git a/test/test_refinement.py b/test/test_refinement.py index 2bc29e89266f4e7a5ae824335d719cb259412ea6..94a4fe0879533163953ae831cb4541b68ca2a951 100644 --- a/test/test_refinement.py +++ b/test/test_refinement.py @@ -27,12 +27,13 @@ from functools import partial import pytest import pyopencl as cl -import pyopencl.clmath # noqa import numpy as np from pyopencl.tools import ( # noqa pytest_generate_tests_for_pyopencl as pytest_generate_tests) +from meshmode.array_context import PyOpenCLArrayContext +from meshmode.dof_array import thaw, flat_norm from meshmode.mesh.generation import ( # noqa generate_icosahedron, generate_box_mesh, make_curve_mesh, ellipse) from meshmode.mesh.refinement.utils import check_nodal_adj_against_geometry @@ -191,6 +192,7 @@ def test_refinement_connection( cl_ctx = ctx_getter() queue = cl.CommandQueue(cl_ctx) + actx = PyOpenCLArrayContext(queue) from meshmode.discretization import Discretization from meshmode.discretization.connection import ( @@ -235,27 +237,27 @@ def test_refinement_connection( factor = 9 for iaxis in range(len(x)): - result = result * cl.clmath.sin(factor * (x[iaxis]/mesh_ext[iaxis])) + result = result * actx.np.sin(factor * (x[iaxis]/mesh_ext[iaxis])) return result - discr = Discretization(cl_ctx, mesh, group_factory(order)) + discr = Discretization(actx, mesh, group_factory(order)) refiner = refiner_cls(mesh) flags = refine_flags(mesh) refiner.refine(flags) connection = make_refinement_connection( - refiner, discr, group_factory(order)) - check_connection(connection) + actx, refiner, discr, group_factory(order)) + check_connection(actx, connection) fine_discr = connection.to_discr - x = discr.nodes().with_queue(queue) - x_fine = fine_discr.nodes().with_queue(queue) + x = thaw(actx, discr.nodes()) + x_fine = thaw(actx, fine_discr.nodes()) f_coarse = f(x) - f_interp = connection(queue, f_coarse).with_queue(queue) - f_true = f(x_fine).with_queue(queue) + f_interp = connection(f_coarse) + f_true = f(x_fine) if visualize == "dots": import matplotlib.pyplot as plt @@ -279,8 +281,7 @@ def test_refinement_connection( ("f_true", f_true), ]) - import numpy.linalg as la - err = la.norm((f_interp - f_true).get(queue), np.inf) + err = flat_norm(f_interp - f_true, np.inf) eoc_rec.add_data_point(h, err) order_slack = 0.5