diff --git a/doc/conf.py b/doc/conf.py index 0cdfb5038..9c6a4522b 100644 --- a/doc/conf.py +++ b/doc/conf.py @@ -12,8 +12,8 @@ # All configuration values have a default; values that are commented out # serve to show the default. -import sys -import os +import sys # noqa: F401 +import os # noqa: F401 # If extensions (or modules to document with autodoc) are in another directory, # add these directories to sys.path here. If the directory is relative to the @@ -23,34 +23,34 @@ # -- General configuration ------------------------------------------------ # If your documentation needs a minimal Sphinx version, state it here. -#needs_sphinx = '1.0' +#needs_sphinx = "1.0" # Add any Sphinx extension module names here, as strings. They can be -# extensions coming with Sphinx (named 'sphinx.ext.*') or your custom +# extensions coming with Sphinx (named "sphinx.ext.*") or your custom # ones. extensions = [ - 'sphinx.ext.autodoc', - 'sphinx.ext.doctest', - 'sphinx.ext.intersphinx', - 'sphinx.ext.mathjax', - 'sphinx.ext.graphviz', + "sphinx.ext.autodoc", + "sphinx.ext.doctest", + "sphinx.ext.intersphinx", + "sphinx.ext.mathjax", + "sphinx.ext.graphviz", ] # Add any paths that contain templates here, relative to this directory. -templates_path = ['_templates'] +templates_path = ["_templates"] # The suffix of source filenames. -source_suffix = '.rst' +source_suffix = ".rst" # The encoding of source files. -#source_encoding = 'utf-8-sig' +#source_encoding = "utf-8-sig" # The master toctree document. -master_doc = 'index' +master_doc = "index" # General information about the project. -project = u'meshmode' -copyright = u'2014, Andreas Klöckner' +project = u"meshmode" +copyright = u"2014, Andreas Klöckner" # The version info for the project you're documenting, acts as replacement for # |version| and |release|, also used in various other places throughout the @@ -58,7 +58,10 @@ # # The short X.Y version. ver_dic = {} -exec(compile(open("../meshmode/version.py").read(), "../meshmode/version.py", 'exec'), ver_dic) +exec( + compile( + open("../meshmode/version.py").read(), "../meshmode/version.py", "exec"), + ver_dic) version = ".".join(str(x) for x in ver_dic["VERSION"]) # The full version, including alpha/beta/rc tags. release = ver_dic["VERSION_TEXT"] @@ -75,7 +78,7 @@ # List of patterns, relative to source directory, that match files and # directories to ignore when looking for source files. -exclude_patterns = ['_build'] +exclude_patterns = ["_build"] # The reST default role (used for this markup: `text`) to use for all # documents. @@ -93,7 +96,7 @@ #show_authors = False # The name of the Pygments (syntax highlighting) style to use. -pygments_style = 'sphinx' +pygments_style = "sphinx" # A list of ignored prefixes for module index sorting. #modindex_common_prefix = [] @@ -114,11 +117,11 @@ } html_sidebars = { - '**': [ - 'about.html', - 'navigation.html', - 'relations.html', - 'searchbox.html', + "**": [ + "about.html", + "navigation.html", + "relations.html", + "searchbox.html", ] } @@ -194,28 +197,28 @@ #html_file_suffix = None # Output file base name for HTML help builder. -htmlhelp_basename = 'meshmodedoc' +htmlhelp_basename = "meshmodedoc" # -- Options for LaTeX output --------------------------------------------- latex_elements = { -# The paper size ('letterpaper' or 'a4paper'). -#'papersize': 'letterpaper', + # The paper size ("letterpaper" or "a4paper"). + #"papersize": "letterpaper", -# The font size ('10pt', '11pt' or '12pt'). -#'pointsize': '10pt', + # The font size ("10pt", "11pt" or "12pt"). + #"pointsize": "10pt", -# Additional stuff for the LaTeX preamble. -#'preamble': '', + # Additional stuff for the LaTeX preamble. + #"preamble": '', } # Grouping the document tree into LaTeX files. List of tuples # (source start file, target name, title, # author, documentclass [howto, manual, or own class]). latex_documents = [ - ('index', 'meshmode.tex', u'meshmode Documentation', - u'Andreas Klöckner', 'manual'), + ("index", "meshmode.tex", u"meshmode Documentation", + u"Andreas Klöckner", "manual"), ] # The name of an image file (relative to this directory) to place at the top of @@ -244,8 +247,8 @@ # One entry per manual page. List of tuples # (source start file, name, description, authors, manual section). man_pages = [ - ('index', 'meshmode', u'meshmode Documentation', - [u'Andreas Klöckner'], 1) + ("index", "meshmode", u"meshmode Documentation", + [u"Andreas Klöckner"], 1) ] # If true, show URL addresses after external links. @@ -258,9 +261,9 @@ # (source start file, target name, title, author, # dir menu entry, description, category) texinfo_documents = [ - ('index', 'meshmode', u'meshmode Documentation', - u'Andreas Klöckner', 'meshmode', 'One line description of project.', - 'Miscellaneous'), + ("index", "meshmode", u"meshmode Documentation", + u"Andreas Klöckner", "meshmode", "One line description of project.", + "Miscellaneous"), ] # Documents to append as an appendix to all manuals. @@ -269,24 +272,25 @@ # If false, no module index is generated. #texinfo_domain_indices = True -# How to display URL addresses: 'footnote', 'no', or 'inline'. -#texinfo_show_urls = 'footnote' +# How to display URL addresses: "footnote", "no", or "inline". +#texinfo_show_urls = "footnote" # If true, do not generate a @detailmenu in the "Top" node's menu. #texinfo_no_detailmenu = False intersphinx_mapping = { - 'https://docs.python.org/': None, - 'https://docs.scipy.org/doc/numpy/': None, - 'https://documen.tician.de/pyopencl': None, - 'https://documen.tician.de/meshpy': None, - 'https://documen.tician.de/modepy': None, - 'https://documen.tician.de/loopy': None, - 'https://documen.tician.de/gmsh_interop': None, - 'https://firedrakeproject.org/': None, - 'https://tisaac.gitlab.io/recursivenodes/': None, - 'https://fenics.readthedocs.io/projects/fiat/en/latest/': None, - 'https://finat.github.io/FInAT/': None, + "https://docs.python.org/3/": None, + "https://numpy.org/doc/stable/": None, + "https://documen.tician.de/pytools": None, + "https://documen.tician.de/pyopencl": None, + "https://documen.tician.de/meshpy": None, + "https://documen.tician.de/modepy": None, + "https://documen.tician.de/loopy": None, + "https://documen.tician.de/gmsh_interop": None, + "https://firedrakeproject.org/": None, + "https://tisaac.gitlab.io/recursivenodes/": None, + "https://fenics.readthedocs.io/projects/fiat/en/latest/": None, + "https://finat.github.io/FInAT/": None, } autoclass_content = "class" diff --git a/examples/simple-dg.py b/examples/simple-dg.py index 6c198d2d8..324f2a925 100644 --- a/examples/simple-dg.py +++ b/examples/simple-dg.py @@ -30,7 +30,7 @@ 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.dof_array import freeze, thaw from meshmode.array_context import PyOpenCLArrayContext, make_loopy_program @@ -146,9 +146,7 @@ def get_connection(self, src, tgt): raise ValueError(f"locations '{src}'->'{tgt}' not understood") def interp(self, src, tgt, vec): - if (isinstance(vec, np.ndarray) - and vec.dtype.char == "O" - and not isinstance(vec, DOFArray)): + if isinstance(vec, np.ndarray): return obj_array_vectorize( lambda el: self.interp(src, tgt, el), vec) @@ -237,9 +235,7 @@ def get_inverse_mass_matrix(self, grp, dtype): return actx.freeze(actx.from_numpy(matrix)) def inverse_mass(self, vec): - if (isinstance(vec, np.ndarray) - and vec.dtype.char == "O" - and not isinstance(vec, DOFArray)): + if isinstance(vec, np.ndarray): return obj_array_vectorize( lambda el: self.inverse_mass(el), vec) @@ -291,11 +287,8 @@ def get_local_face_mass_matrix(self, afgrp, volgrp, dtype): return actx.freeze(actx.from_numpy(matrix)) def face_mass(self, vec): - if (isinstance(vec, np.ndarray) - and vec.dtype.char == "O" - and not isinstance(vec, DOFArray)): - return obj_array_vectorize( - lambda el: self.face_mass(el), vec) + if isinstance(vec, np.ndarray): + return obj_array_vectorize(lambda el: self.face_mass(el), vec) @memoize_in(self, "face_mass_knl") def knl(): diff --git a/meshmode/array_context.py b/meshmode/array_context.py index 72d48201f..f9a490037 100644 --- a/meshmode/array_context.py +++ b/meshmode/array_context.py @@ -25,7 +25,6 @@ import loopy as lp from loopy.version import MOST_RECENT_LANGUAGE_VERSION from pytools import memoize_method -from pytools.obj_array import obj_array_vectorized_n_args __doc__ = """ .. autofunction:: make_loopy_program @@ -58,7 +57,7 @@ def __init__(self, array_context): self._array_context = array_context def __getattr__(self, name): - def f(*args): + def loopy_implemented_elwise_func(*args): actx = self._array_context # FIXME: Maybe involve loopy type inference? result = actx.empty(args[0].shape, args[0].dtype) @@ -68,14 +67,15 @@ def f(*args): **{"inp%d" % i: arg for i, arg in enumerate(args)}) return result - return obj_array_vectorized_n_args(f) + from meshmode.dof_array import obj_or_dof_array_vectorized_n_args + return obj_or_dof_array_vectorized_n_args(loopy_implemented_elwise_func) - @obj_array_vectorized_n_args def conjugate(self, x): # NOTE: conjugate distribute over object arrays, but it looks for a # `conjugate` ufunc, while some implementations only have the shorter # `conj` (e.g. cl.array.Array), so this should work for everybody - return x.conj() + from meshmode.dof_array import obj_or_dof_array_vectorize + return obj_or_dof_array_vectorize(lambda obj: obj.conj(), x) conj = conjugate @@ -225,20 +225,21 @@ def thaw(self, array): # {{{ PyOpenCLArrayContext class _PyOpenCLFakeNumpyNamespace(_BaseFakeNumpyNamespace): - @obj_array_vectorized_n_args def maximum(self, x, y): import pyopencl.array as cl_array - return cl_array.maximum(x, y) + from meshmode.dof_array import obj_or_dof_array_vectorize_n_args + return obj_or_dof_array_vectorize_n_args(cl_array.maximum, x, y) - @obj_array_vectorized_n_args def minimum(self, x, y): import pyopencl.array as cl_array - return cl_array.minimum(x, y) + from meshmode.dof_array import obj_or_dof_array_vectorize_n_args + return obj_or_dof_array_vectorize_n_args(cl_array.minimum, x, y) - @obj_array_vectorized_n_args def where(self, criterion, then, else_): import pyopencl.array as cl_array - return cl_array.if_positive(criterion != 0, then, else_) + from meshmode.dof_array import obj_or_dof_array_vectorize_n_args + return obj_or_dof_array_vectorize_n_args( + cl_array.if_positive, criterion != 0, then, else_) class PyOpenCLArrayContext(ArrayContext): diff --git a/meshmode/discretization/__init__.py b/meshmode/discretization/__init__.py index d1ec5a6b5..f164d5d72 100644 --- a/meshmode/discretization/__init__.py +++ b/meshmode/discretization/__init__.py @@ -237,9 +237,9 @@ def _new_array(self, actx, creation_func, dtype=None): else: dtype = np.dtype(dtype) - return _DOFArray.from_list(actx, [ + return _DOFArray(actx, tuple( creation_func(shape=(grp.nelements, grp.nunit_dofs), dtype=dtype) - for grp in self.groups]) + for grp in self.groups)) def empty(self, actx: ArrayContext, dtype=None): """Return an empty :class:`~meshmode.dof_array.DOFArray`. @@ -297,11 +297,11 @@ def get_mat(grp): return mat - return _DOFArray.from_list(actx, [ + return _DOFArray(actx, tuple( actx.call_loopy( prg(), diff_mat=actx.from_numpy(get_mat(grp)), vec=vec[grp.index] )["result"] - for grp in self.groups]) + for grp in self.groups)) @memoize_method def quad_weights(self): @@ -316,14 +316,14 @@ def prg(): "result[iel,idof] = weights[idof]", name="quad_weights") - return _DOFArray.from_list(None, [ + return _DOFArray(None, tuple( actx.freeze( actx.call_loopy( prg(), weights=actx.from_numpy(grp.weights), nelements=grp.nelements, )["result"]) - for grp in self.groups]) + for grp in self.groups)) @memoize_method def nodes(self): @@ -348,7 +348,7 @@ def prg(): name="nodes") return make_obj_array([ - _DOFArray.from_list(None, [ + _DOFArray(None, tuple( actx.freeze( actx.call_loopy( prg(), @@ -356,7 +356,7 @@ def prg(): grp.from_mesh_interp_matrix()), nodes=actx.from_numpy(grp.mesh_el_group.nodes[iaxis]) )["result"]) - for grp in self.groups]) + for grp in self.groups)) for iaxis in range(self.ambient_dim)]) # vim: fdm=marker diff --git a/meshmode/discretization/visualization.py b/meshmode/discretization/visualization.py index 8eb339475..2426b25cd 100644 --- a/meshmode/discretization/visualization.py +++ b/meshmode/discretization/visualization.py @@ -22,7 +22,7 @@ import numpy as np from pytools import memoize_method, Record -from meshmode.dof_array import DOFArray, flatten, thaw +from meshmode.dof_array import flatten, thaw __doc__ = """ @@ -73,9 +73,7 @@ def separate_by_real_and_imag(names_and_fields, real_only): def resample_to_numpy(conn, vec): - if (isinstance(vec, np.ndarray) - and vec.dtype.char == "O" - and not isinstance(vec, DOFArray)): + if isinstance(vec, np.ndarray): from pytools.obj_array import obj_array_vectorize return obj_array_vectorize(lambda x: resample_to_numpy(conn, x), vec) diff --git a/meshmode/dof_array.py b/meshmode/dof_array.py index c1bde0163..6cff7bbbd 100644 --- a/meshmode/dof_array.py +++ b/meshmode/dof_array.py @@ -22,17 +22,26 @@ import operator import numpy as np -from typing import Optional, Iterable, Any +from typing import Optional, Iterable, Any, Tuple, Union from functools import partial +from numbers import Number +import operator as op +import decorator from pytools import single_valued, memoize_in -from pytools.obj_array import obj_array_vectorize, obj_array_vectorize_n_args +from pytools.obj_array import obj_array_vectorize from meshmode.array_context import ArrayContext, make_loopy_program __doc__ = """ .. autoclass:: DOFArray + +.. autofunction:: obj_or_dof_array_vectorize +.. autofunction:: obj_or_dof_array_vectorized +.. autofunction:: obj_or_dof_array_vectorize_n_args +.. autofunction:: obj_or_dof_array_vectorized_n_args + .. autofunction:: thaw .. autofunction:: freeze .. autofunction:: flatten @@ -43,9 +52,8 @@ # {{{ 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 DOFArray: + """This array type holds degree-of-freedom arrays for use with :class:`~meshmode.discretization.Discretization`, with one entry in the :class:`DOFArray` for each :class:`~meshmode.discretization.ElementGroupBase`. @@ -56,9 +64,7 @@ class DOFArray(np.ndarray): of the associated group. ``ndofs_per_element`` is typically, but not necessarily, the same as :attr:`~meshmode.discretization.ElementGroupBase.nunit_dofs` - of the associated group. - This array is derived from :class:`numpy.ndarray` with dtype object ("an - object array"). The entries in this array are further arrays managed by + of the associated group. The entries in this array are further arrays managed by :attr:`array_context`. One main purpose of this class is to describe the data structure, @@ -75,7 +81,10 @@ class DOFArray(np.ndarray): The (assumed uniform) :class:`numpy.dtype` of the group arrays contained in this instance. - .. automethod:: from_list + .. automethod:: __len__ + .. automethod:: __getitem__ + + This object supports arithmetic, comparisons, and logic operators. .. note:: @@ -83,28 +92,23 @@ class DOFArray(np.ndarray): ``<=``, ``>=``. (:mod:`numpy` object arrays containing arrays do not.) """ - # Follows https://numpy.org/devdocs/user/basics.subclassing.html - - def __new__(cls, actx: Optional[ArrayContext], input_array): + def __init__(self, actx: Optional[ArrayContext], data: Tuple[Any]): 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") + if not isinstance(data, tuple): + raise TypeError("'data' argument must be a tuple") - result.array_context = actx - return result + self.array_context = actx + self._data = data - def __array_finalize__(self, obj): - if obj is None: - return - self.array_context = getattr(obj, "array_context", None) + # Tell numpy that we would like to do our own array math, thank you very much. + # (numpy arrays have priority 0.) + __array_priority__ = 10 @property def entry_dtype(self): - return single_valued(subary.dtype for subary in self.flat) + return single_valued(subary.dtype for subary in self._data) @classmethod def from_list(cls, actx: Optional[ArrayContext], res_list) -> "DOFArray": @@ -114,96 +118,183 @@ def from_list(cls, actx: Optional[ArrayContext], res_list) -> "DOFArray": :arg actx: If *None*, the arrays in *res_list* must be :meth:`~meshmode.array_context.ArrayContext.thaw`\ ed. """ + from warnings import warn + warn("DOFArray.from_list is deprecated and will disappear in 2021.", + DeprecationWarning, stacklevel=2) if not (actx is None or isinstance(actx, ArrayContext)): raise TypeError("actx must be of type ArrayContext") - ngroups = len(res_list) + return cls(actx, tuple(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] + # {{{ sequence protocol - return result + def __len__(self): + return len(self._data) - # {{{ work around numpy failing to compare obj arrays of arrays + def __getitem__(self, i): + return self._data[i] - def _comparison(self, operator_func, other): - from numbers import Number - if isinstance(other, DOFArray): - return obj_array_vectorize_n_args(operator_func, self, other) + def __iter__(self): + return iter(self._data) - elif isinstance(other, Number): - return obj_array_vectorize( - lambda self_entry: operator_func(self_entry, other), - self) + # }}} + @property + def shape(self): + return (len(self),) + + def _like_me(self, data): + return DOFArray(self.array_context, tuple(data)) + + def _bop(self, f, op1, op2): + """Broadcasting logic for a generic binary operator.""" + if isinstance(op1, DOFArray) and isinstance(op2, DOFArray): + if len(op1._data) != len(op2._data): + raise ValueError("DOFArray objects in binary operator must have " + f"same length, got {len(op1._data)} and {len(op2._data)}") + return self._like_me([ + f(op1_i, op2_i) + for op1_i, op2_i in zip(op1._data, op2._data)]) + elif isinstance(op1, DOFArray) and isinstance(op2, Number): + return self._like_me([f(op1_i, op2) for op1_i in op1._data]) + elif isinstance(op1, Number) and isinstance(op2, DOFArray): + return self._like_me([f(op1, op2_i) for op2_i in op2._data]) else: - # fall back to "best effort" (i.e. likley failure) - return operator_func(self, other) + return NotImplemented + + def __add__(self, arg): return self._bop(op.add, self, arg) # noqa: E704 + def __radd__(self, arg): return self._bop(op.add, arg, self) # noqa: E704 + def __sub__(self, arg): return self._bop(op.sub, self, arg) # noqa: E704 + def __rsub__(self, arg): return self._bop(op.sub, arg, self) # noqa: E704 + def __mul__(self, arg): return self._bop(op.mul, self, arg) # noqa: E704 + def __rmul__(self, arg): return self._bop(op.mul, arg, self) # noqa: E704 + def __truediv__(self, arg): return self._bop(op.truediv, self, arg) # noqa: E704 + def __rtruediv__(self, arg): return self._bop(op.truediv, arg, self) # noqa: E704, E501 + def __pow__(self, arg): return self._bop(op.pow, self, arg) # noqa: E704 + def __rpow__(self, arg): return self._bop(op.pow, arg, self) # noqa: E704 + def __mod__(self, arg): return self._bop(op.mod, self, arg) # noqa: E704 + def __rmod__(self, arg): return self._bop(op.mod, arg, self) # noqa: E704 + def __divmod__(self, arg): return self._bop(divmod, self, arg) # noqa: E704 + def __rdivmod__(self, arg): return self._bop(divmod, arg, self) # noqa: E704 + + def __pos__(self): return self # noqa: E704 + def __neg__(self): return self._like_me([-self_i for self_i in self._data]) # noqa: E704, E501 + def __abs__(self): return self._like_me([abs(self_i) for self_i in self._data]) # noqa: E704, E501 + + def conj(self): return self._like_me([self_i.conj() for self_i in self._data]) # noqa: E704, E501 + @property + def real(self): return self._like_me([self_i.real for self_i in self._data]) # noqa: E704, E501 + @property + def imag(self): return self._like_me([self_i.imag for self_i in self._data]) # noqa: E704, E501 - def __lt__(self, other): - return self._comparison(operator.lt, other) + def __eq__(self, arg): return self._bop(op.eq, self, arg) # noqa: E704 + def __ne__(self, arg): return self._bop(op.ne, self, arg) # noqa: E704 + def __lt__(self, arg): return self._bop(op.lt, self, arg) # noqa: E704 + def __gt__(self, arg): return self._bop(op.gt, self, arg) # noqa: E704 + def __le__(self, arg): return self._bop(op.le, self, arg) # noqa: E704 + def __ge__(self, arg): return self._bop(op.ge, self, arg) # noqa: E704 - def __gt__(self, other): - return self._comparison(operator.gt, other) + def __and__(self, arg): return self._bop(operator.and_, self, arg) # noqa: E704 + def __xor__(self, arg): return self._bop(operator.xor, self, arg) # noqa: E704 + def __or__(self, arg): return self._bop(operator.or_, self, arg) # noqa: E704 + def __rand__(self, arg): return self._bop(operator.and_, arg, self) # noqa: E704 + def __rxor__(self, arg): return self._bop(operator.xor, arg, self) # noqa: E704 + def __ror__(self, arg): return self._bop(operator.or_, arg, self) # noqa: E704 - def __le__(self, other): - return self._comparison(operator.le, other) + # bit shifts unimplemented for now - def __ge__(self, other): - return self._comparison(operator.ge, other) +# }}} - # }}} -# }}} +def obj_or_dof_array_vectorize(f, ary): + r""" + Works like :func:`~pytools.obj_array.obj_array_vectorize`, but also + for :class:`DOFArray`\ s. + """ + + if isinstance(ary, DOFArray): + return ary._like_me([f(ary_i) for ary_i in ary._data]) + else: + return obj_array_vectorize(f, ary) + + +obj_or_dof_array_vectorized = decorator.decorator(obj_or_dof_array_vectorize) + + +def obj_or_dof_array_vectorize_n_args(f, *args): + r"""Apply the function *f* elementwise to all entries of any + object arrays or :class:`DOFArray`\ s in *args*. All such arrays are expected + to have the same shape (but this is not checked). + Equivalent to an appropriately-looped execution of:: + + result[idx] = f(obj_array_arg1[idx], arg2, obj_array_arg3[idx]) + + Return an array of the same shape as the arguments consisting of the + return values of *f*. + + Works like :func:`~pytools.obj_array.obj_array_vectorize_n_args`, but also + for :class:`DOFArray`\ s. + """ + dofarray_arg_indices = [ + i for i, arg in enumerate(args) + if isinstance(arg, DOFArray)] + + if not dofarray_arg_indices: + from pytools.obj_array import obj_array_vectorize_n_args + return obj_array_vectorize_n_args(f, *args) + + leading_da_index = dofarray_arg_indices[0] + + template_ary = args[leading_da_index] + result = [] + new_args = list(args) + for igrp in range(len(template_ary)): + for arg_i in dofarray_arg_indices: + new_args[arg_i] = args[arg_i][igrp] + result.append(f(*new_args)) + + return DOFArray(template_ary.array_context, tuple(result)) + + +obj_or_dof_array_vectorized_n_args = decorator.decorator( + obj_or_dof_array_vectorize_n_args) -def thaw(actx: ArrayContext, ary: np.ndarray) -> np.ndarray: +def thaw(actx: ArrayContext, ary: Union[DOFArray, 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)): + if isinstance(ary, np.ndarray): 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 - ]) + return DOFArray(actx, tuple(actx.thaw(subary) for subary in ary)) -def freeze(ary: np.ndarray) -> np.ndarray: +def freeze(ary: Union[DOFArray, 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)): + if isinstance(ary, np.ndarray): 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 - ]) + return DOFArray(None, tuple( + ary.array_context.freeze(subary) for subary in ary)) -def flatten(ary: np.ndarray) -> Any: +def flatten(ary: Union[DOFArray, np.ndarray]) -> Any: 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`. @@ -214,9 +305,7 @@ def flatten(ary: np.ndarray) -> Any: Vectorizes over object arrays of :class:`DOFArray`\ s. """ - if (isinstance(ary, np.ndarray) - and ary.dtype.char == "O" - and not isinstance(ary, DOFArray)): + if isinstance(ary, np.ndarray): return obj_array_vectorize(flatten, ary) group_sizes = [grp_ary.shape[0] * grp_ary.shape[1] for grp_ary in ary] @@ -240,15 +329,13 @@ def prg(): return result -def unflatten(actx: ArrayContext, discr, ary, +def unflatten(actx: ArrayContext, discr, ary: Union[Any, np.ndarray], ndofs_per_element_per_group: Optional[Iterable[int]] = None) -> 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)): + if isinstance(ary, np.ndarray): return obj_array_vectorize( lambda subary: unflatten( actx, discr, subary, ndofs_per_element_per_group), @@ -276,7 +363,7 @@ def prg(): group_starts = np.cumsum([0] + group_sizes) - return DOFArray.from_list(actx, [ + return DOFArray(actx, tuple( actx.call_loopy( prg(), grp_start=grp_start, ary=ary, @@ -286,7 +373,7 @@ def prg(): for grp_start, grp, ndofs_per_element in zip( group_starts, discr.groups, - ndofs_per_element_per_group)]) + ndofs_per_element_per_group))) def flat_norm(ary: DOFArray, ord=2): diff --git a/meshmode/interop/nodal_dg.py b/meshmode/interop/nodal_dg.py index 82e5e1468..98fc278e4 100644 --- a/meshmode/interop/nodal_dg.py +++ b/meshmode/interop/nodal_dg.py @@ -156,7 +156,7 @@ def pull_dof_array( ) -> meshmode.dof_array.DOFArray: ary = self.octave.pull(name).T - return meshmode.dof_array.DOFArray.from_list(actx, [actx.from_numpy(ary)]) + return meshmode.dof_array.DOFArray(actx, (actx.from_numpy(ary),)) def download_nodal_dg_if_not_present(path="nodal-dg"): diff --git a/test/test_array.py b/test/test_array.py new file mode 100644 index 000000000..07da804a5 --- /dev/null +++ b/test/test_array.py @@ -0,0 +1,176 @@ +__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 + +import meshmode # noqa: F401 +from meshmode.array_context import ( # noqa + pytest_generate_tests_for_pyopencl_array_context + as pytest_generate_tests) + +from meshmode.discretization import Discretization +from meshmode.discretization.poly_element import PolynomialWarpAndBlendGroupFactory +from meshmode.dof_array import flatten, unflatten + +import logging +logger = logging.getLogger(__name__) + + +def test_array_context_np_workalike(actx_factory): + actx = actx_factory() + + from meshmode.mesh.generation import generate_regular_rect_mesh + mesh = generate_regular_rect_mesh( + a=(-0.5,)*2, b=(0.5,)*2, n=(8,)*2, order=3) + + discr = Discretization(actx, mesh, PolynomialWarpAndBlendGroupFactory(3)) + + for sym_name, n_args in [ + ("sin", 1), + ("exp", 1), + ("arctan2", 2), + ("minimum", 2), + ("maximum", 2), + ("where", 3), + ("conj", 1), + ]: + args = [np.random.randn(discr.ndofs) for i in range(n_args)] + ref_result = getattr(np, sym_name)(*args) + + actx_args = [unflatten(actx, discr, actx.from_numpy(arg)) for arg in args] + + actx_result = actx.to_numpy( + flatten(getattr(actx.np, sym_name)(*actx_args))) + + assert np.allclose(actx_result, ref_result) + + +def test_dof_array_arithmetic_same_as_numpy(actx_factory): + actx = actx_factory() + + from meshmode.mesh.generation import generate_regular_rect_mesh + mesh = generate_regular_rect_mesh( + a=(-0.5,)*2, b=(0.5,)*2, n=(3,)*2, order=1) + + discr = Discretization(actx, mesh, PolynomialWarpAndBlendGroupFactory(3)) + + import operator + from pytools import generate_nonnegative_integer_tuples_below as gnitb + from random import uniform, randrange + for op_func, n_args, use_integers in [ + (operator.add, 2, False), + (operator.sub, 2, False), + (operator.mul, 2, False), + (operator.truediv, 2, False), + (operator.pow, 2, False), + # FIXME pyopencl.Array doesn't do mod. + #(operator.mod, 2, True), + #(operator.mod, 2, False), + # FIXME: Two outputs + #(divmod, 2, False), + + (operator.and_, 2, True), + (operator.xor, 2, True), + (operator.or_, 2, True), + + (operator.le, 2, False), + (operator.ge, 2, False), + (operator.lt, 2, False), + (operator.gt, 2, False), + (operator.eq, 2, True), + (operator.ne, 2, True), + + (operator.pos, 1, False), + (operator.neg, 1, False), + (operator.abs, 1, False), + + (lambda ary: ary.real, 1, False), + (lambda ary: ary.imag, 1, False), + ]: + for is_array_flags in gnitb(2, n_args): + if sum(is_array_flags) == 0: + # all scalars, no need to test + continue + args = [ + (0.5+np.random.rand(discr.ndofs) + if not use_integers else + np.random.randint(3, 200, discr.ndofs)) + + if is_array_flag else + (uniform(0.5, 2) + if not use_integers + else randrange(3, 200)) + for is_array_flag in is_array_flags] + ref_result = op_func(*args) + + actx_args = [ + unflatten(actx, discr, actx.from_numpy(arg)) + if isinstance(arg, np.ndarray) else arg + for arg in args] + + actx_result = actx.to_numpy(flatten(op_func(*actx_args))) + + assert np.allclose(actx_result, ref_result) + + +def test_dof_array_comparison(actx_factory): + actx = actx_factory() + + from meshmode.mesh.generation import generate_regular_rect_mesh + mesh = generate_regular_rect_mesh( + a=(-0.5,)*2, b=(0.5,)*2, n=(8,)*2, order=3) + + discr = Discretization(actx, mesh, PolynomialWarpAndBlendGroupFactory(3)) + + import operator + for op in [ + operator.lt, + operator.le, + operator.gt, + operator.ge, + ]: + np_arg = np.random.randn(discr.ndofs) + arg = unflatten(actx, discr, actx.from_numpy(np_arg)) + zeros = discr.zeros(actx) + + comp = op(arg, zeros) + np_comp = actx.to_numpy(flatten(comp)) + assert np.array_equal(np_comp, op(np_arg, 0)) + + comp = op(arg, 0) + np_comp = actx.to_numpy(flatten(comp)) + assert np.array_equal(np_comp, op(np_arg, 0)) + + comp = op(0, arg) + np_comp = actx.to_numpy(flatten(comp)) + assert np.array_equal(np_comp, op(0, np_arg)) + + +if __name__ == "__main__": + import sys + if len(sys.argv) > 1: + exec(sys.argv[1]) + else: + from pytest import main + main([__file__]) + +# vim: fdm=marker diff --git a/test/test_firedrake_interop.py b/test/test_firedrake_interop.py index f8189f526..045978c89 100644 --- a/test/test_firedrake_interop.py +++ b/test/test_firedrake_interop.py @@ -666,7 +666,7 @@ def test_to_fd_idempotency(ctx_factory, mm_mesh, fspace_degree): mm_unique = discr.zeros(actx, dtype=dtype) unique_vals = np.arange(np.size(mm_unique[0]), dtype=dtype) mm_unique[0].set(unique_vals.reshape(mm_unique[0].shape)) - mm_unique_copy = DOFArray.from_list(actx, [mm_unique[0].copy()]) + mm_unique_copy = DOFArray(actx, (mm_unique[0].copy(),)) # Test for idempotency mm->fd->mm fdrake_unique = fdrake_connection.from_meshmode(mm_unique) diff --git a/test/test_meshmode.py b/test/test_meshmode.py index cf90befb7..f64513f7e 100644 --- a/test/test_meshmode.py +++ b/test/test_meshmode.py @@ -40,7 +40,7 @@ LegendreGaussLobattoTensorProductGroupFactory ) from meshmode.mesh import Mesh, BTAG_ALL -from meshmode.dof_array import thaw, flat_norm, flatten, unflatten +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 @@ -191,7 +191,7 @@ def test_visualizers(actx_factory, dim, group_factory): discr = Discretization(actx, mesh, discr_group_factory(target_order)) nodes = thaw(actx, discr.nodes()) - f = actx.np.sqrt(sum(nodes**2)) + f = actx.np.sqrt(sum(nodes**2)) + 1j*nodes[0] from meshmode.discretization.visualization import make_visualizer vis = make_visualizer(actx, discr, target_order) @@ -1589,74 +1589,6 @@ def test_mesh_multiple_groups(actx_factory, ambient_dim, visualize=False): discr.num_reference_derivative(ref_axes, x[0]) -def test_array_context_np_workalike(actx_factory): - actx = actx_factory() - - from meshmode.mesh.generation import generate_regular_rect_mesh - mesh = generate_regular_rect_mesh( - a=(-0.5,)*2, b=(0.5,)*2, n=(8,)*2, order=3) - - from meshmode.discretization import Discretization - from meshmode.discretization.poly_element import \ - PolynomialWarpAndBlendGroupFactory as GroupFactory - discr = Discretization(actx, mesh, GroupFactory(3)) - - for sym_name, n_args in [ - ("sin", 1), - ("exp", 1), - ("arctan2", 2), - ("minimum", 2), - ("maximum", 2), - ("where", 3), - ("conj", 1), - ]: - args = [np.random.randn(discr.ndofs) for i in range(n_args)] - ref_result = getattr(np, sym_name)(*args) - - actx_args = [unflatten(actx, discr, actx.from_numpy(arg)) for arg in args] - - actx_result = actx.to_numpy( - flatten(getattr(actx.np, sym_name)(*actx_args))) - - assert np.allclose(actx_result, ref_result) - - -def test_dof_array_comparison(actx_factory): - actx = actx_factory() - - from meshmode.mesh.generation import generate_regular_rect_mesh - mesh = generate_regular_rect_mesh( - a=(-0.5,)*2, b=(0.5,)*2, n=(8,)*2, order=3) - - from meshmode.discretization import Discretization - from meshmode.discretization.poly_element import \ - PolynomialWarpAndBlendGroupFactory as GroupFactory - discr = Discretization(actx, mesh, GroupFactory(3)) - - import operator - for op in [ - operator.lt, - operator.le, - operator.gt, - operator.ge, - ]: - np_arg = np.random.randn(discr.ndofs) - arg = unflatten(actx, discr, actx.from_numpy(np_arg)) - zeros = discr.zeros(actx) - - comp = op(arg, zeros) - np_comp = actx.to_numpy(flatten(comp)) - assert np.array_equal(np_comp, op(np_arg, 0)) - - comp = op(arg, 0) - np_comp = actx.to_numpy(flatten(comp)) - assert np.array_equal(np_comp, op(np_arg, 0)) - - comp = op(0, arg) - np_comp = actx.to_numpy(flatten(comp)) - assert np.array_equal(np_comp, op(0, np_arg)) - - if __name__ == "__main__": import sys if len(sys.argv) > 1: