From 5bfbf7c50938dbad8399e497f8c86e2cc7327055 Mon Sep 17 00:00:00 2001 From: Andreas Kloeckner Date: Thu, 8 Oct 2020 17:36:49 -0500 Subject: [PATCH 1/8] Stub out interface for nodal-DG interop --- doc/interop.rst | 4 ++ meshmode/interop/nodal_dg.py | 88 ++++++++++++++++++++++++++++++++++++ 2 files changed, 92 insertions(+) create mode 100644 meshmode/interop/nodal_dg.py diff --git a/doc/interop.rst b/doc/interop.rst index 7f6acb60..66537f00 100644 --- a/doc/interop.rst +++ b/doc/interop.rst @@ -4,6 +4,10 @@ Interoperability with Other Discretization Packages Functionality in this subpackage helps import and export data to/from other pieces of software, typically PDE solvers. +Nodal DG +-------- + +.. automodule:: meshmode.interop.nodal_dg Firedrake --------- diff --git a/meshmode/interop/nodal_dg.py b/meshmode/interop/nodal_dg.py new file mode 100644 index 00000000..47fcd223 --- /dev/null +++ b/meshmode/interop/nodal_dg.py @@ -0,0 +1,88 @@ +"""Provides interoperability with the `Matlab/Octave Codes +`__ complementing the +book "Nodal Discontinuous Galerkin Methods" by Jan Hesthaven +and Tim Warburton (Springer, 2008). + +.. autoclass:: NodalDGContext +""" + +__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 meshmode.mesh +import meshmode.discretization +import meshmode.dof_array +import meshmode.array_context + + +class NodalDGContext(object): + """Should be used as a context manager to ensure proper cleanup. + + .. automethod:: __init__ + .. automethod:: set_mesh + .. automethod:: get_discr + .. automethod:: push_dof_array + .. automethod:: pull_dof_array + """ + + def __init__(self, path): + """ + :arg path: The path to the ``Codes1.1`` folder of the nodal DG codes. + """ + self.path = path + self.octave = None + + def __enter__(self): + import oct2py + self.octave = oct2py.Oct2Py() + import os.path + self.octave.feval(os.path.join(self.path, "mypath.m")) + + def __exit__(self, exc_type, exc_val, exc_tb): + self.octave.kill_octave() + + def set_mesh(self, mesh: meshmode.mesh.Mesh): + """Set the mesh information in the nodal DG Octave instance to + the one given by *mesh*. + + High-order information is silently ignored. + """ + pass + + def get_discr(self) -> meshmode.discretization.Discretization: + """Get a discretization with nodes exactly matching the ones used + by the book code. + + The returned discretization contains a new :class:`~meshmode.mesh.Mesh` + object constructed from the global Octave state. + """ + + def push_dof_array(self, name, ary: meshmode.dof_array.DOFArray): + """ + """ + pass + + def pull_dof_array( + self, actx: meshmode.array_context.ArrayContext, name + ) -> meshmode.dof_array.DOFArray: + pass -- GitLab From 905f1ca3c08ad30597fc388ed92ba317e8624675 Mon Sep 17 00:00:00 2001 From: Andreas Kloeckner Date: Fri, 9 Oct 2020 00:44:30 -0500 Subject: [PATCH 2/8] Finish implementing, testing nodal-dg interop --- meshmode/discretization/poly_element.py | 42 ++++++++++++ meshmode/interop/nodal_dg.py | 89 ++++++++++++++++++++++--- 2 files changed, 121 insertions(+), 10 deletions(-) diff --git a/meshmode/discretization/poly_element.py b/meshmode/discretization/poly_element.py index 3651a670..cf257613 100644 --- a/meshmode/discretization/poly_element.py +++ b/meshmode/discretization/poly_element.py @@ -40,6 +40,7 @@ Group types .. autoclass:: PolynomialWarpAndBlendElementGroup .. autoclass:: PolynomialRecursiveNodesElementGroup .. autoclass:: PolynomialEquidistantSimplexElementGroup +.. autoclass:: PolynomialGivenNodesElementGroup .. autoclass:: LegendreGaussLobattoTensorProductElementGroup .. autoclass:: EquidistantTensorProductElementGroup @@ -52,6 +53,7 @@ Group factories .. autoclass:: PolynomialWarpAndBlendGroupFactory .. autoclass:: PolynomialRecursiveNodesGroupFactory .. autoclass:: PolynomialEquidistantSimplexGroupFactory +.. autoclass:: PolynomialGivenNodesGroupFactory .. autoclass:: LegendreGaussLobattoTensorProductGroupFactory """ @@ -312,6 +314,33 @@ class PolynomialEquidistantSimplexElementGroup(_MassMatrixQuadratureElementGroup assert dim2 == dim return result + +class PolynomialGivenNodesElementGroup(_MassMatrixQuadratureElementGroup): + """Elemental discretization with a number of nodes matching the number of + polynomials in :math:`P^k`, hence usable for differentiation and + interpolation. Uses nodes given by the user. + """ + def __init__(self, mesh_el_group, order, unit_nodes, index): + super().__init__(mesh_el_group, order, index) + self._unit_nodes = unit_nodes + + @property + def unit_nodes(self): + dim2, nunit_nodes = self._unit_nodes.shape + + if dim2 != self.mesh_el_group.dim: + raise ValueError("unit nodes supplied to " + "PolynomialGivenNodesElementGroup do not have expected " + "dimensionality") + + if nunit_nodes != len(self.basis()): + raise ValueError("unit nodes supplied to " + "PolynomialGivenNodesElementGroup do not have expected " + "node count for provided order") + + return self._unit_nodes + + # }}} @@ -441,6 +470,19 @@ class PolynomialEquidistantSimplexGroupFactory(HomogeneousOrderBasedGroupFactory mesh_group_class = _MeshSimplexElementGroup group_class = PolynomialEquidistantSimplexElementGroup + +class PolynomialGivenNodesGroupFactory(HomogeneousOrderBasedGroupFactory): + def __init__(self, order, unit_nodes): + self.order = order + self.unit_nodes = unit_nodes + + def __call__(self, mesh_el_group, index): + if not isinstance(mesh_el_group, _MeshSimplexElementGroup): + raise TypeError("only mesh element groups of type '%s' " + "are supported" % _MeshSimplexElementGroup.__name__) + + return PolynomialGivenNodesElementGroup( + mesh_el_group, self.order, self.unit_nodes, index) # }}} diff --git a/meshmode/interop/nodal_dg.py b/meshmode/interop/nodal_dg.py index 47fcd223..a97969bd 100644 --- a/meshmode/interop/nodal_dg.py +++ b/meshmode/interop/nodal_dg.py @@ -29,6 +29,7 @@ THE SOFTWARE. """ +import numpy as np import meshmode.mesh import meshmode.discretization import meshmode.dof_array @@ -55,34 +56,102 @@ class NodalDGContext(object): def __enter__(self): import oct2py self.octave = oct2py.Oct2Py() - import os.path - self.octave.feval(os.path.join(self.path, "mypath.m")) + self.octave.eval(f'cd "{self.path}"') + self.octave.eval("mypath") + return self def __exit__(self, exc_type, exc_val, exc_tb): - self.octave.kill_octave() + self.octave.exit() - def set_mesh(self, mesh: meshmode.mesh.Mesh): + REF_AXES = ["r", "s", "t"] + AXES = ["x", "y", "z"] + + def set_mesh(self, mesh: meshmode.mesh.Mesh, order): """Set the mesh information in the nodal DG Octave instance to the one given by *mesh*. - High-order information is silently ignored. + The mesh must only have a single element group of simplices. + + .. warning:: + + High-order geometryinformation is currently silently ignored. """ - pass + if len(mesh.groups) != 1: + raise ValueError("mesh must have exactly one element group") + + elgrp, = mesh.groups + + self.octave.eval(f"Globals{mesh.dim}D;") + self.octave.push("Nv", mesh.nvertices) + self.octave.push("K", mesh.nelements) + for ax in range(mesh.ambient_dim): + self.octave.push(f"V{self.AXES[ax].upper()}", mesh.vertices[ax]) - def get_discr(self) -> meshmode.discretization.Discretization: + self.octave.push(f"V{self.AXES[ax].upper()}", mesh.vertices[ax]) + self.octave.push("EToV", elgrp.vertex_indices+1) + + self.octave.push("N", order) + + self.octave.eval(f"StartUp{mesh.dim}D;") + + def get_discr(self, actx) -> meshmode.discretization.Discretization: """Get a discretization with nodes exactly matching the ones used - by the book code. + by the nodal-DG code. The returned discretization contains a new :class:`~meshmode.mesh.Mesh` object constructed from the global Octave state. """ + # find dim as number of vertices in the simplex - 1 + etov_size = self.octave.eval("size(EToV)", verbose=False) + dim = int(etov_size[0, 1]-1) + + if dim == 1: + unit_nodes = self.octave.eval(f"JacobiGL(0, 0, N)", verbose=False).T + else: + unit_nodes_arrays = self.octave.eval( + f"Nodes{dim}D(N)", nout=dim, verbose=False) + + equilat_to_biunit_func_name = ( + "".join(self.AXES[:dim] + ["to"] + self.REF_AXES[:dim])) + + unit_nodes_arrays = self.octave.feval( + equilat_to_biunit_func_name, *unit_nodes_arrays, + nout=dim, verbose=False) + + unit_nodes = np.array([a.reshape(-1) for a in unit_nodes_arrays]) + + vertices = np.array([ + self.octave.pull(f"V{self.AXES[ax].upper()}").reshape(-1) + for ax in range(dim)]) + nodes = np.array([self.octave.pull(self.AXES[ax]).T for ax in range(dim)]) + vertex_indices = (self.octave.pull("EToV")).astype(np.int32)-1 + + from meshmode.mesh import Mesh, SimplexElementGroup + order = int(self.octave.pull("N")) + egroup = SimplexElementGroup( + order, + vertex_indices=vertex_indices, + nodes=nodes, + unit_nodes=unit_nodes) + + mesh = Mesh(vertices=vertices, groups=[egroup], is_conforming=True) + + from meshmode.discretization import Discretization + from meshmode.discretization.poly_element import ( + PolynomialGivenNodesGroupFactory) + return Discretization(actx, mesh, + PolynomialGivenNodesGroupFactory(order, unit_nodes)) def push_dof_array(self, name, ary: meshmode.dof_array.DOFArray): """ """ - pass + grp_array, = ary + ary = ary.array_context.to_numpy(grp_array) + self.octave.push(name, ary.T) def pull_dof_array( self, actx: meshmode.array_context.ArrayContext, name ) -> meshmode.dof_array.DOFArray: - pass + ary = self.octave.pull(name).T + + return meshmode.dof_array.DOFArray.from_list(actx, [actx.from_numpy(ary)]) -- GitLab From ca38f4830cb7c5091688034f9924781ef07c356c Mon Sep 17 00:00:00 2001 From: Andreas Kloeckner Date: Fri, 9 Oct 2020 00:45:53 -0500 Subject: [PATCH 3/8] Placate flake8 --- meshmode/interop/nodal_dg.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/meshmode/interop/nodal_dg.py b/meshmode/interop/nodal_dg.py index a97969bd..8eb72154 100644 --- a/meshmode/interop/nodal_dg.py +++ b/meshmode/interop/nodal_dg.py @@ -106,7 +106,7 @@ class NodalDGContext(object): dim = int(etov_size[0, 1]-1) if dim == 1: - unit_nodes = self.octave.eval(f"JacobiGL(0, 0, N)", verbose=False).T + unit_nodes = self.octave.eval("JacobiGL(0, 0, N)", verbose=False).T else: unit_nodes_arrays = self.octave.eval( f"Nodes{dim}D(N)", nout=dim, verbose=False) -- GitLab From af04f2e6155fbd30167aa14ffde13eaf86630cae Mon Sep 17 00:00:00 2001 From: Andreas Kloeckner Date: Fri, 9 Oct 2020 00:47:35 -0500 Subject: [PATCH 4/8] Add nodal-dg interop test file --- test/test_interop.py | 81 ++++++++++++++++++++++++++++++++++++++++++++ 1 file changed, 81 insertions(+) create mode 100644 test/test_interop.py diff --git a/test/test_interop.py b/test/test_interop.py new file mode 100644 index 00000000..bbbf46d4 --- /dev/null +++ b/test/test_interop.py @@ -0,0 +1,81 @@ +__copyright__ = "Copyright (C) 2014 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 pytest +import pyopencl as cl + +from meshmode.array_context import ( # noqa + pytest_generate_tests_for_pyopencl_array_context + as pytest_generate_tests, + PyOpenCLArrayContext) +from meshmode.dof_array import flat_norm, thaw + +import logging +logger = logging.getLogger(__name__) + + +def acf(): + context = cl._csc() + queue = cl.CommandQueue(context) + return PyOpenCLArrayContext(queue) + + +@pytest.mark.parametrize("dim", [1, 2, 3]) +def test_nodal_dg_interop(actx_factory, dim): + actx = actx_factory() + + order = 4 + + from meshmode.mesh.generation import generate_regular_rect_mesh + mesh = generate_regular_rect_mesh( + a=(-0.5,)*dim, b=(0.5,)*dim, n=(8,)*dim, order=order) + + from meshmode.interop.nodal_dg import NodalDGContext + with NodalDGContext("./nodal-dg/Codes1.1") as ndgctx: + ndgctx.set_mesh(mesh, order=order) + + discr = ndgctx.get_discr(actx) + + for ax in range(dim): + x_ax = ndgctx.pull_dof_array(actx, ndgctx.AXES[ax]) + err = flat_norm(x_ax-discr.nodes()[ax], np.inf) + assert err < 1e-15 + + n0 = thaw(actx, discr.nodes()[0]) + + ndgctx.push_dof_array("n0", n0) + n0_2 = ndgctx.pull_dof_array(actx, "n0") + + assert flat_norm(n0 - n0_2, np.inf) < 1e-15 + + +if __name__ == "__main__": + import sys + if len(sys.argv) > 1: + exec(sys.argv[1]) + else: + from pytest import main + main([__file__]) + +# vim: fdm=marker -- GitLab From 00d4d7ee7f55aad9c80e4023e783d32ff9871f20 Mon Sep 17 00:00:00 2001 From: Andreas Kloeckner Date: Fri, 9 Oct 2020 17:01:59 -0500 Subject: [PATCH 5/8] CI tweaks to enable testing nodal-DG interop --- .github/workflows/ci.yml | 1 + .gitlab-ci.yml | 15 ++++++++------- .test-conda-env-py3.yml | 3 +++ 3 files changed, 12 insertions(+), 7 deletions(-) diff --git a/.github/workflows/ci.yml b/.github/workflows/ci.yml index 7ccf2775..7f8eeb6a 100644 --- a/.github/workflows/ci.yml +++ b/.github/workflows/ci.yml @@ -30,6 +30,7 @@ jobs: - name: "Main Script" run: | CONDA_ENVIRONMENT=.test-conda-env-py3.yml + (cd test; git clone https://github.com/tcew/nodal-dg) curl -L -O -k https://gitlab.tiker.net/inducer/ci-support/raw/master/build-and-test-py-project-within-miniconda.sh . ./build-and-test-py-project-within-miniconda.sh diff --git a/.gitlab-ci.yml b/.gitlab-ci.yml index 3c29f131..80b52c39 100644 --- a/.gitlab-ci.yml +++ b/.gitlab-ci.yml @@ -1,11 +1,12 @@ Python 3 POCL: - script: - - export PY_EXE=python3 - - export PYOPENCL_TEST=portable:pthread - # cython is here because pytential (for now, for TS) depends on it - - export EXTRA_INSTALL="pybind11 cython numpy mako mpi4py" - - curl -L -O -k https://gitlab.tiker.net/inducer/ci-support/raw/master/build-and-test-py-project.sh - - ". ./build-and-test-py-project.sh" + script: | + export PY_EXE=python3 + export PYOPENCL_TEST=portable:pthread + (cd test; git clone https://github.com/tcew/nodal-dg) + # cython is here because pytential (for now, for TS) depends on it + export EXTRA_INSTALL="pybind11 cython numpy mako mpi4py oct2py" + curl -L -O -k https://gitlab.tiker.net/inducer/ci-support/raw/master/build-and-test-py-project.sh + . ./build-and-test-py-project.sh tags: - python3 - pocl diff --git a/.test-conda-env-py3.yml b/.test-conda-env-py3.yml index 05d64de7..5da5fe79 100644 --- a/.test-conda-env-py3.yml +++ b/.test-conda-env-py3.yml @@ -13,6 +13,9 @@ dependencies: - islpy - gmsh +# Needed for nodal-DG interop +- oct2py + # for Pytential - cython -- GitLab From 3164f68741b450ae5ef4c4e1ebfcfd9ed2dd5de0 Mon Sep 17 00:00:00 2001 From: Andreas Kloeckner Date: Sat, 10 Oct 2020 23:04:56 -0500 Subject: [PATCH 6/8] Install octave in Conda to allow testing nodal-dg interop --- .test-conda-env-py3.yml | 1 + 1 file changed, 1 insertion(+) diff --git a/.test-conda-env-py3.yml b/.test-conda-env-py3.yml index 5da5fe79..0b15e85b 100644 --- a/.test-conda-env-py3.yml +++ b/.test-conda-env-py3.yml @@ -15,6 +15,7 @@ dependencies: # Needed for nodal-DG interop - oct2py +- octave # for Pytential - cython -- GitLab From ba566d2509fa030f373561c2205487f0fc18f41c Mon Sep 17 00:00:00 2001 From: Andreas Kloeckner Date: Sat, 10 Oct 2020 23:09:09 -0500 Subject: [PATCH 7/8] Also clone nodal DG for Gitlab Conda CI --- .gitlab-ci.yml | 1 + 1 file changed, 1 insertion(+) diff --git a/.gitlab-ci.yml b/.gitlab-ci.yml index 80b52c39..e4968466 100644 --- a/.gitlab-ci.yml +++ b/.gitlab-ci.yml @@ -77,6 +77,7 @@ Python 3 POCL Firedrake Examples: Python 3 Conda: script: | CONDA_ENVIRONMENT=.test-conda-env-py3.yml + (cd test; git clone https://github.com/tcew/nodal-dg) curl -L -O -k https://gitlab.tiker.net/inducer/ci-support/raw/master/build-and-test-py-project-within-miniconda.sh . ./build-and-test-py-project-within-miniconda.sh tags: -- GitLab From aac0ae71f3bcb9064e3da624ace079dd31d2d01d Mon Sep 17 00:00:00 2001 From: Andreas Kloeckner Date: Sat, 10 Oct 2020 23:44:47 -0500 Subject: [PATCH 8/8] Github CI: Install Octave via apt, not via conda-forge --- .github/workflows/ci.yml | 3 +++ .test-conda-env-py3.yml | 1 - 2 files changed, 3 insertions(+), 1 deletion(-) diff --git a/.github/workflows/ci.yml b/.github/workflows/ci.yml index 7f8eeb6a..79fab479 100644 --- a/.github/workflows/ci.yml +++ b/.github/workflows/ci.yml @@ -29,6 +29,9 @@ jobs: - uses: actions/checkout@v2 - name: "Main Script" run: | + sudo apt update + sudo apt install octave + CONDA_ENVIRONMENT=.test-conda-env-py3.yml (cd test; git clone https://github.com/tcew/nodal-dg) curl -L -O -k https://gitlab.tiker.net/inducer/ci-support/raw/master/build-and-test-py-project-within-miniconda.sh diff --git a/.test-conda-env-py3.yml b/.test-conda-env-py3.yml index 0b15e85b..5da5fe79 100644 --- a/.test-conda-env-py3.yml +++ b/.test-conda-env-py3.yml @@ -15,7 +15,6 @@ dependencies: # Needed for nodal-DG interop - oct2py -- octave # for Pytential - cython -- GitLab