From 0179a830a98b554003ca2bd3ed72997db7a71ff9 Mon Sep 17 00:00:00 2001 From: Alexandru Fikl Date: Sat, 13 Jul 2024 15:27:57 +0300 Subject: [PATCH 1/5] use Mesh.copy instead of replace --- meshmode/mesh/processing.py | 3 +-- test/test_mesh.py | 3 +-- test/test_meshmode.py | 2 +- 3 files changed, 3 insertions(+), 5 deletions(-) diff --git a/meshmode/mesh/processing.py b/meshmode/mesh/processing.py index 1c23fbc0..9e03f253 100644 --- a/meshmode/mesh/processing.py +++ b/meshmode/mesh/processing.py @@ -1646,8 +1646,7 @@ def not_none(vi: Optional[np.ndarray]) -> np.ndarray: new_vertex_indices = np.cumsum(used_flags, dtype=mesh.vertex_id_dtype) - 1 new_vertex_indices[~used_flags] = -1 - return replace( - mesh, + return mesh.copy( vertices=mesh.vertices[:, used_flags], groups=tuple( replace(grp, vertex_indices=new_vertex_indices[grp.vertex_indices]) diff --git a/test/test_mesh.py b/test/test_mesh.py index 6b0d4c87..c471e095 100644 --- a/test/test_mesh.py +++ b/test/test_mesh.py @@ -267,8 +267,7 @@ def test_remove_unused_vertices(): assert mesh.vertices is not None - mesh2 = replace( - mesh, + mesh2 = mesh.copy( vertices=np.concatenate([np.zeros((3, 1)), mesh.vertices], axis=1), groups=tuple( replace( diff --git a/test/test_meshmode.py b/test/test_meshmode.py index 526e7a2e..f264b646 100644 --- a/test/test_meshmode.py +++ b/test/test_meshmode.py @@ -889,7 +889,7 @@ def test_mesh_without_vertices(actx_factory): groups = [ replace(grp, nodes=grp.nodes, vertex_indices=None) for grp in mesh.groups] - mesh = make_mesh(None, groups, is_conforming=False) + mesh = make_mesh(None, groups, is_conforming=None) # try refining it from meshmode.mesh.refinement import refine_uniformly From 6cef90d8372e6429feb52ca4898bf06671db1bcc Mon Sep 17 00:00:00 2001 From: Alexandru Fikl Date: Sat, 13 Jul 2024 15:35:59 +0300 Subject: [PATCH 2/5] remove last call to actx.np.norm on a DOFArray --- test/test_array.py | 23 +++++++++++------------ 1 file changed, 11 insertions(+), 12 deletions(-) diff --git a/test/test_array.py b/test/test_array.py index f6a3b209..15af9e6d 100644 --- a/test/test_array.py +++ b/test/test_array.py @@ -28,6 +28,7 @@ from arraycontext import ( dataclass_array_container, + flatten, pytest_generate_tests_for_array_contexts, with_container_arithmetic, ) @@ -142,33 +143,31 @@ def _get_test_containers(actx, ambient_dim=2): @pytest.mark.parametrize("ord", [2, np.inf]) def test_container_norm(actx_factory, ord): actx = actx_factory() - c_test = _get_test_containers(actx) - # {{{ actx.np.linalg.norm + # {{{ flat_norm from pytools.obj_array import make_obj_array c = MyContainer(name="hey", mass=1, momentum=make_obj_array([2, 3]), enthalpy=5) c_obj_ary = make_obj_array([c, c]) - n1 = actx.np.linalg.norm(c_obj_ary, ord) + n1 = flat_norm(c_obj_ary, ord) n2 = np.linalg.norm([1, 2, 3, 5]*2, ord) assert abs(n1 - n2) < 1e-12 - # }}} - - # {{{ flat_norm - - # check nested vs actx.np.linalg.norm - assert actx.to_numpy(abs( - flat_norm(c_test[1], ord=ord) - - actx.np.linalg.norm(c_test[1], ord=ord))) < 1e-12 - # check nested container with only Numbers (and no actx) assert abs(flat_norm(c_obj_ary, ord=ord) - n2) < 1.0e-12 assert abs( flat_norm(np.array([1, 1], dtype=object), ord=ord) - np.linalg.norm([1, 1], ord=ord)) < 1.0e-12 + + # check nested + n1 = actx.to_numpy(flat_norm(c_test[1], ord=ord)) + n2 = np.linalg.norm([ + np.linalg.norm(actx.to_numpy(flatten(ary, actx)), ord=ord) for ary in c_test[1] + ], ord=ord) + assert abs(n1 - n2) < 1e-12 + # }}} # }}} From 752a5f7e1655ea13f68edd8182154cc110edcb44 Mon Sep 17 00:00:00 2001 From: Alexandru Fikl Date: Sat, 13 Jul 2024 15:57:54 +0300 Subject: [PATCH 3/5] port inverse_mass_matrix to take a Basis --- examples/simple-dg.py | 4 +--- meshmode/discretization/poly_element.py | 14 +++++++------- 2 files changed, 8 insertions(+), 10 deletions(-) diff --git a/examples/simple-dg.py b/examples/simple-dg.py index 138b7d65..f8752698 100644 --- a/examples/simple-dg.py +++ b/examples/simple-dg.py @@ -237,9 +237,7 @@ def face_jacobian(self, where): @memoize_method def get_inverse_mass_matrix(self, grp, dtype): import modepy as mp - matrix = mp.inverse_mass_matrix( - grp.basis_obj().functions, - grp.unit_nodes) + matrix = mp.inverse_mass_matrix(grp.basis_obj(), grp.unit_nodes) actx = self._setup_actx return actx.freeze(actx.from_numpy(matrix)) diff --git a/meshmode/discretization/poly_element.py b/meshmode/discretization/poly_element.py index da6ec2b8..6fcfae5f 100644 --- a/meshmode/discretization/poly_element.py +++ b/meshmode/discretization/poly_element.py @@ -126,7 +126,7 @@ def mass_matrix(grp: InterpolatoryElementGroupBase) -> np.ndarray: assert grp.is_orthonormal_basis() return mp.mass_matrix( - grp.basis_obj().functions, + grp.basis_obj(), grp.unit_nodes) @@ -288,11 +288,11 @@ def quadrature_rule(self): class _MassMatrixQuadratureElementGroup(PolynomialSimplexElementGroupBase): @memoize_method def quadrature_rule(self): - basis_fcts = self.basis_obj().functions + basis = self.basis_obj() nodes = self._interp_nodes - mass_matrix = mp.mass_matrix(basis_fcts, nodes) + mass_matrix = mp.mass_matrix(basis, nodes) weights = np.dot(mass_matrix, - np.ones(len(basis_fcts))) + np.ones(len(basis.functions))) return mp.Quadrature(nodes, weights, exact_to=self.order) @property @@ -571,11 +571,11 @@ def basis_obj(self): @memoize_method def quadrature_rule(self): - basis_fcts = self._basis.functions + basis = self._basis nodes = self._nodes - mass_matrix = mp.mass_matrix(basis_fcts, nodes) + mass_matrix = mp.mass_matrix(basis, nodes) weights = np.dot(mass_matrix, - np.ones(len(basis_fcts))) + np.ones(len(basis.functions))) return mp.Quadrature(nodes, weights, exact_to=self.order) @property From d610a61b9709b3ed6387783e72c19abd17b57d76 Mon Sep 17 00:00:00 2001 From: Alexandru Fikl Date: Sat, 13 Jul 2024 16:04:24 +0300 Subject: [PATCH 4/5] fix some deprecations in examples --- examples/moving-geometry.py | 2 +- examples/plot-connectivity.py | 2 +- examples/simple-dg.py | 4 +++- 3 files changed, 5 insertions(+), 3 deletions(-) diff --git a/examples/moving-geometry.py b/examples/moving-geometry.py index efec9fe3..81fe086d 100644 --- a/examples/moving-geometry.py +++ b/examples/moving-geometry.py @@ -245,7 +245,7 @@ def source(t, x): if __name__ == "__main__": cl_ctx = cl.create_some_context() queue = cl.CommandQueue(cl_ctx) - actx = PyOpenCLArrayContext(queue) + actx = PyOpenCLArrayContext(queue, force_device_scalars=True) from pytools import ProcessTimer for _ in range(1): diff --git a/examples/plot-connectivity.py b/examples/plot-connectivity.py index ce27d5db..7d3bbd59 100644 --- a/examples/plot-connectivity.py +++ b/examples/plot-connectivity.py @@ -11,7 +11,7 @@ def main(): cl_ctx = cl.create_some_context() queue = cl.CommandQueue(cl_ctx) - actx = PyOpenCLArrayContext(queue) + actx = PyOpenCLArrayContext(queue, force_device_scalars=True) from meshmode.mesh.generation import ( # noqa: F401 generate_icosahedron, diff --git a/examples/simple-dg.py b/examples/simple-dg.py index f8752698..52ae4547 100644 --- a/examples/simple-dg.py +++ b/examples/simple-dg.py @@ -439,7 +439,9 @@ def bump(actx, discr, t=0): / source_width**2)) -@with_container_arithmetic(bcast_obj_array=True, rel_comparison=True) +@with_container_arithmetic(bcast_obj_array=True, + rel_comparison=True, + _cls_has_array_context_attr=True) @dataclass_array_container @dataclass(frozen=True) class WaveState: From 3ea05d6e8a0e19f70242e6e3153aa48bd3c249ba Mon Sep 17 00:00:00 2001 From: Alexandru Fikl Date: Sat, 13 Jul 2024 16:41:37 +0300 Subject: [PATCH 5/5] poly_element: return _interp_nodes directly for MassMatrix groups --- meshmode/discretization/poly_element.py | 5 +++++ 1 file changed, 5 insertions(+) diff --git a/meshmode/discretization/poly_element.py b/meshmode/discretization/poly_element.py index 6fcfae5f..9d511e74 100644 --- a/meshmode/discretization/poly_element.py +++ b/meshmode/discretization/poly_element.py @@ -295,6 +295,11 @@ def quadrature_rule(self): np.ones(len(basis.functions))) return mp.Quadrature(nodes, weights, exact_to=self.order) + @property + @memoize_method + def unit_nodes(self): + return self._interp_nodes + @property @abstractmethod def _interp_nodes(self):