Tensor product 1d nodes & basis functions #384

a-alveyblanc · 2023-08-04T18:40:36Z

Without access to the 1D nodes and 1D basis functions used in the tensor product, we are not able to exploit the tensor product structure to reduce the cost required for operator evaluation, e.g. applying a differentiation operator to function data.

This PR:

Adds functionality to TensorProductElementGroupBase to store and retrieve 1D basis functions and 1D nodes used in the tensor product
Updates subclasses accordingly

Depends on:

TensorProductBasis: drop orth_weight argument modepy#72

alexfikl

Left a few comments that will hopefully simplify things, but otherwise looks good!

alexfikl · 2023-08-05T17:43:33Z

meshmode/discretization/poly_element.py

+    def __init__(self, mesh_el_group, order, index=None, *, basis, basis_1d,
+                 unit_nodes, unit_nodes_1d):


It feels a bit weird to pass in both (basis, unit_nodes) and (basis_1d, unit_nodes_1d) here since we can construct the first one from the second. Maybe move the construction of the nd ones here?

Fixed this for 1D nodes in my latest push. modepy doesn't seem to have a function like mp.tensor_product_nodes for basis functions, so still working on that.

I don't know if one is needed. You should be able to just construct a TensorProductBasis like here (for a more general case, you'd just have a tensor product of mesh_el_group.dim 1D functions)
https://github.com/inducer/modepy/blob/0fa7355ebe979a4706a7da2ae15b7517d229b8b3/modepy/modes.py#L1219-L1232

Would that work?

meshmode/discretization/poly_element.py

alexfikl · 2023-08-05T17:44:28Z

meshmode/discretization/poly_element.py

@@ -549,6 +559,10 @@ def quadrature_rule(self):
                         np.ones(len(basis_fcts)))
        return mp.Quadrature(nodes, weights, exact_to=self.order)

+    @property
+    def unit_nodes_1d(self):
+        return np.array([self._nodes_1d])


Why is the np.array needed?

Added a comment to clarify this. In short, modepy expects a nested array as an argument when evaluating the basis functions.

Maybe do this via reshape instead (to clarify that it's already an array).

inducer · 2023-08-07T14:47:06Z

.pylintrc-local.yml

This wants to stay.

inducer · 2023-08-07T14:54:10Z

meshmode/discretization/poly_element.py

        """
        :arg basis: a :class:`modepy.TensorProductBasis`.
-        :arg unit_nodes: unit nodes for the tensor product, obtained by
-            using :func:`modepy.tensor_product_nodes`, for example.
+        :arg basis_1d: a representation of the 1D basis from which the tensor


Let's make the 1D bases accessible from non-underscored attributes in modepy.

inducer · 2023-08-07T14:56:08Z

meshmode/discretization/poly_element.py

@@ -549,6 +559,10 @@ def quadrature_rule(self):
                         np.ones(len(basis_fcts)))
        return mp.Quadrature(nodes, weights, exact_to=self.order)

+    @property
+    def unit_nodes_1d(self):
+        return np.array([self._nodes_1d])


Maybe do this via reshape instead (to clarify that it's already an array).

inducer · 2023-08-07T14:57:18Z

meshmode/discretization/poly_element.py

+        unit_nodes = mp.tensor_product_nodes(
+                [self.unit_nodes_1d[0]] * self.mesh_el_group.dim)
+
+        if unit_nodes.shape[0] != self.mesh_el_group.dim:


Maybe call self.unit_nodes on construction so that this error arises at that point (and not arbitrarily later)?

a-alveyblanc · 2023-09-06T16:31:33Z

Converted to draft since there's a bit more work I'd like to do other than making 1D nodes accessible.

modepy will have some updated functionality for returning 1D operators in the case of tensor products
meshmode will need some updating to reflect this

LMK if you have any questions

inducer · 2023-11-13T17:04:26Z

meshmode/discretization/poly_element.py

+        # FIXME?? there are cases where `basis` is a _SimplexONB object, not a
+        # TensorProductBasis object. in those cases, it seems self._bases_1d is
+        # unused. leaving as-is for now.
+        if isinstance(basis, mp.TensorProductBasis):
+            self._bases_1d = basis.bases[0]


Instead: If it's a 1D affair (check on both basis and nodes before making that decision), and the basis is not TP, wrap the basis in a TensorProductBasis for consistency.

inducer · 2023-11-13T17:05:09Z

meshmode/discretization/poly_element.py

+    def __init__(self, mesh_el_group, order, index=None, *, basis,
+                 unit_nodes):


Try type-annotating?

a-alveyblanc · 2023-11-14T14:20:52Z

This is ready for a final look I think. Only thing I'm unsure about is the type annotations I added to TensorProductElementGroupBase

inducer

LGTM with those two fixes.

meshmode/discretization/poly_element.py

inducer · 2023-11-14T21:11:48Z

I think the Firedrake bits are failing because https://pypi.org/project/pocl-binary-distribution/ (aka https://github.com/isuruf/ocl-wheels/) does not have an LLVM new enough to understand EPYC CPUs.

inducer · 2023-11-14T21:13:02Z

I tried to do a rebuild, but the container pull is failing with some sort of AWS signature message that I don't know how to fix. Maybe https://quay.io is broken ATM?

a-alveyblanc · 2023-11-14T21:17:10Z

Other tests are passing on AMD Epyc, though. Is there a large difference between other tests and Firedrake? Grudge tests passed on Epyc: test/test_dt_utils.py::test_wave_dt_estimate[<PytatoPyOpenCLArrayContext for <pyopencl.Device 'cpu-AMD EPYC 7763 64-Core Processor' on 'Portable Computing Language'>>-2-1] (grabbed some random warning output from the Grudge tests)

inducer · 2023-11-14T21:30:39Z

They are! It depends on how pocl is installed. By and large, we install pocl from conda, and that works fine. For Firedrake, we install pocl via pip (yes, that's possible... Isuru worked some magic). That's the bit that's broken.

alexfikl · 2023-11-14T21:34:23Z

They are! It depends on how pocl is installed. By and large, we install pocl from conda, and that works fine. For Firedrake, we install pocl via pip (yes, that's possible... Isuru worked some magic). That's the bit that's broken.

Doesn't it install it from the Ubuntu repos?

meshmode/.ci/install-for-firedrake.sh

Line 10 in 9e942e6

sudo apt install -y pocl-opencl-icd ocl-icd-opencl-dev

inducer · 2023-11-15T17:11:31Z

Good catch! Thanks for spotting that. That makes for an obvious thing to try: Use pocl from pypi. :)

inducer · 2023-11-15T17:19:37Z

Would you look at that. ;) Firedrake did run on EPYC and succeeded.

…he tensor product Add functionality to store and retrieve 1D basis functions of the tensor product Return 1d nodes as a nested array to match shape expected by modepy Construct tensor product nodes from 1D nodes in TensorProductElementGroupBase Remove files created by running pylint Fix a small typo Minor change Get 1D and ND basis by passing a callable to TPGroupBase, probably not final Change formatting on docstring Mesh TP element groups: claim to be non-affine by default Change docstring Take advantage of changes in modepy to grab 1D components of TP Add FIXME regarding basis not being a TensorProductBasisObject Minor changes Minor documentation change Remove trailing whitespace Update handling of non-TensorProductBasis objects in TensorProductElementGroup Refactor Minor change Add type annotations Stop supplying orth weight to TP basis ValueError -> NotImplementedError Minor fixes Fully type TensorProductElementGroupBase constructor

a-alveyblanc mentioned this pull request Aug 4, 2023

Tensor product differentiation operators inducer/grudge#313

Closed

1 task

alexfikl reviewed Aug 5, 2023

View reviewed changes

inducer reviewed Aug 7, 2023

View reviewed changes

a-alveyblanc marked this pull request as draft September 6, 2023 16:30

inducer reviewed Nov 13, 2023

View reviewed changes

a-alveyblanc marked this pull request as ready for review November 14, 2023 14:20

inducer approved these changes Nov 14, 2023

View reviewed changes

meshmode/discretization/poly_element.py Outdated Show resolved Hide resolved

meshmode/discretization/poly_element.py Outdated Show resolved Hide resolved

inducer force-pushed the tensor-product-1d-nodes-and-1d-basis branch from 4442850 to 4dc68ba Compare November 14, 2023 16:23

a-alveyblanc and others added 2 commits November 15, 2023 11:27

Firedrake CI: Install pocl via PyPI

f07a64c

inducer force-pushed the tensor-product-1d-nodes-and-1d-basis branch from 7ec0efc to f07a64c Compare November 15, 2023 17:27

inducer enabled auto-merge (rebase) November 15, 2023 17:28

Install graphviz and pdf2svg via apt for docs

f6a3600

inducer merged commit 052d11e into inducer:main Nov 15, 2023
12 checks passed

Provide feedback

Saved searches

Use saved searches to filter your results more quickly

Tensor product 1d nodes & basis functions #384

Tensor product 1d nodes & basis functions #384

a-alveyblanc commented Aug 4, 2023 •

edited by inducer

Loading

alexfikl left a comment

alexfikl Aug 5, 2023

a-alveyblanc Aug 7, 2023 •

edited

Loading

alexfikl Aug 8, 2023 •

edited

Loading

alexfikl Aug 5, 2023

a-alveyblanc Aug 7, 2023

inducer Aug 7, 2023

inducer Aug 7, 2023

inducer Aug 7, 2023

inducer Aug 7, 2023

inducer Aug 7, 2023

a-alveyblanc commented Sep 6, 2023 •

edited

Loading

inducer Nov 13, 2023

inducer Nov 13, 2023

a-alveyblanc commented Nov 14, 2023 •

edited

Loading

inducer left a comment

inducer commented Nov 14, 2023

inducer commented Nov 14, 2023

a-alveyblanc commented Nov 14, 2023 •

edited

Loading

inducer commented Nov 14, 2023

alexfikl commented Nov 14, 2023

inducer commented Nov 15, 2023

inducer commented Nov 15, 2023

		def __init__(self, mesh_el_group, order, index=None, *, basis, basis_1d,
		unit_nodes, unit_nodes_1d):

		def __init__(self, mesh_el_group, order, index=None, *, basis,
		unit_nodes):

Tensor product 1d nodes & basis functions #384

Tensor product 1d nodes & basis functions #384

Conversation

a-alveyblanc commented Aug 4, 2023 • edited by inducer Loading

alexfikl left a comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

a-alveyblanc Aug 7, 2023 • edited Loading

Choose a reason for hiding this comment

alexfikl Aug 8, 2023 • edited Loading

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

a-alveyblanc commented Sep 6, 2023 • edited Loading

Choose a reason for hiding this comment

Choose a reason for hiding this comment

a-alveyblanc commented Nov 14, 2023 • edited Loading

inducer left a comment

Choose a reason for hiding this comment

inducer commented Nov 14, 2023

inducer commented Nov 14, 2023

a-alveyblanc commented Nov 14, 2023 • edited Loading

inducer commented Nov 14, 2023

alexfikl commented Nov 14, 2023

inducer commented Nov 15, 2023

inducer commented Nov 15, 2023

a-alveyblanc commented Aug 4, 2023 •

edited by inducer

Loading

a-alveyblanc Aug 7, 2023 •

edited

Loading

alexfikl Aug 8, 2023 •

edited

Loading

a-alveyblanc commented Sep 6, 2023 •

edited

Loading

a-alveyblanc commented Nov 14, 2023 •

edited

Loading

a-alveyblanc commented Nov 14, 2023 •

edited

Loading