diff --git a/.github/workflows/ci.yml b/.github/workflows/ci.yml index 7ccf277560b917963d5a1ba57b402d4cc9837c9d..79fab4791d9da2dfc918c29cf2c6e484185275e6 100644 --- a/.github/workflows/ci.yml +++ b/.github/workflows/ci.yml @@ -29,7 +29,11 @@ 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 . ./build-and-test-py-project-within-miniconda.sh diff --git a/.gitlab-ci.yml b/.gitlab-ci.yml index 3c29f131cb3f9bcbe4c01e042ab446770951e9fe..e4968466ec98e13d97e724929256ff9123727333 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 @@ -76,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: diff --git a/.test-conda-env-py3.yml b/.test-conda-env-py3.yml index 05d64de7d6dafb8c7537a9fa660a6458b558727f..5da5fe79d6dc2f10b4bd41ebf7e1a7a01cb14648 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 diff --git a/doc/interop.rst b/doc/interop.rst index 7f6acb60daf9585afa1fc2ca85c599da6a9d781b..66537f00787f817710a849836a6d479715e13063 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/discretization/poly_element.py b/meshmode/discretization/poly_element.py index 3651a6707212aeb908977cd783f62b278502094b..cf2576137a434451986068cfc18415f386a68b0e 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 new file mode 100644 index 0000000000000000000000000000000000000000..8eb721546ff927239e8e2506a962983ed0c2047f --- /dev/null +++ b/meshmode/interop/nodal_dg.py @@ -0,0 +1,157 @@ +"""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 numpy as np +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() + self.octave.eval(f'cd "{self.path}"') + self.octave.eval("mypath") + return self + + def __exit__(self, exc_type, exc_val, exc_tb): + self.octave.exit() + + 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*. + + The mesh must only have a single element group of simplices. + + .. warning:: + + High-order geometryinformation is currently silently ignored. + """ + 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]) + + 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 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("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): + """ + """ + 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: + ary = self.octave.pull(name).T + + return meshmode.dof_array.DOFArray.from_list(actx, [actx.from_numpy(ary)]) diff --git a/test/test_interop.py b/test/test_interop.py new file mode 100644 index 0000000000000000000000000000000000000000..bbbf46d47fe81a88d82cfbbd186d2b03cef7076b --- /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