From 8f8b310e632044e9b7456ed68f6a0a2e6d62089a Mon Sep 17 00:00:00 2001 From: Andreas Kloeckner Date: Fri, 24 May 2024 17:01:43 -0500 Subject: [PATCH] Add Gmsh{2,3}DMeshBuilder Also, refactor test MeshBuilder to not pretend to act as source of truth on mesh order. --- test/mesh_data.py | 96 ++++++++++++++++++++++----------- test/test_dt_utils.py | 31 +++++------ test/test_grudge.py | 115 +++++++++++++++++++++++++--------------- test/test_reductions.py | 29 +++++----- 4 files changed, 163 insertions(+), 108 deletions(-) diff --git a/test/mesh_data.py b/test/mesh_data.py index 0ccc369a1..c71950bb9 100644 --- a/test/mesh_data.py +++ b/test/mesh_data.py @@ -1,35 +1,51 @@ +from abc import abstractmethod, ABC +from typing import ClassVar, Hashable, Optional, Sequence +from meshmode.mesh import Mesh +from meshmode.mesh.io import read_gmsh import numpy as np import meshmode.mesh.generation as mgen -class MeshBuilder: - order = 4 - mesh_order = None +class MeshBuilder(ABC): + resolutions: ClassVar[Sequence[Hashable]] + ambient_dim: ClassVar[int] - def __init__(self, **kwargs): - for k, v in kwargs.items(): - setattr(self, k, v) + @abstractmethod + def get_mesh( + self, + resolution: Hashable, + mesh_order: Optional[int] = None + ) -> Mesh: + ... - if self.mesh_order is None: - self.mesh_order = self.order - @property - def ambient_dim(self): - raise NotImplementedError +class _GmshMeshBuilder(MeshBuilder): + resolutions = [None] - @property - def resolutions(self): - raise NotImplementedError + def __init__(self, filename: str) -> None: + self._mesh_fn = filename + + def get_mesh(self, resolution, mesh_order=None) -> Mesh: + assert resolution is None + assert mesh_order is None + return read_gmsh(self._mesh_fn, force_ambient_dim=self.ambient_dim) + + +class GmshMeshBuilder2D(_GmshMeshBuilder): + ambient_dim = 2 - def get_mesh(self, resolution, mesh_order): - raise NotImplementedError + +class GmshMeshBuilder3D(_GmshMeshBuilder): + ambient_dim = 3 class Curve2DMeshBuilder(MeshBuilder): ambient_dim = 2 resolutions = [16, 32, 64, 128] - def get_mesh(self, resolution, mesh_order): + def get_mesh(self, resolution, mesh_order=None): + if mesh_order is None: + mesh_order = 4 return mgen.make_curve_mesh( self.curve_fn, # pylint: disable=no-member np.linspace(0.0, 1.0, resolution + 1), @@ -37,8 +53,9 @@ def get_mesh(self, resolution, mesh_order): class EllipseMeshBuilder(Curve2DMeshBuilder): - radius = 3.1 - aspect_ratio = 2.0 + def __init__(self, radius=3.1, aspect_ratio=2): + self.radius = radius + self.aspect_ratio = aspect_ratio @property def curve_fn(self): @@ -58,9 +75,13 @@ class SphereMeshBuilder(MeshBuilder): ambient_dim = 3 resolutions = [0, 1, 2, 3] - radius = 1.0 - def get_mesh(self, resolution, mesh_order): + radius: float + + def __init__(self, radius=1): + self.radius = radius + + def get_mesh(self, resolution, mesh_order=4): from meshmode.mesh.generation import generate_sphere return generate_sphere(self.radius, order=mesh_order, uniform_refinement_rounds=resolution) @@ -69,13 +90,16 @@ def get_mesh(self, resolution, mesh_order): class SpheroidMeshBuilder(MeshBuilder): ambient_dim = 3 - mesh_order = 4 resolutions = [0, 1, 2, 3] - radius = 1.0 - aspect_ratio = 2.0 + radius: float + aspect_ratio: float + + def __init__(self, radius=1, aspect_ratio=2): + self.radius = radius + self.aspect_ratio = aspect_ratio - def get_mesh(self, resolution, mesh_order): + def get_mesh(self, resolution, mesh_order=4): from meshmode.mesh.generation import generate_sphere mesh = generate_sphere(self.radius, order=mesh_order, uniform_refinement_rounds=resolution) @@ -84,16 +108,14 @@ def get_mesh(self, resolution, mesh_order): return affine_map(mesh, A=np.diag([1.0, 1.0, self.aspect_ratio])) -class BoxMeshBuilder(MeshBuilder): - ambient_dim = 2 - - mesh_order = 1 +class _BoxMeshBuilderBase(MeshBuilder): resolutions = [4, 8, 16] + mesh_order = 1 a = (-0.5, -0.5, -0.5) b = (+0.5, +0.5, +0.5) - def get_mesh(self, resolution, mesh_order): + def get_mesh(self, resolution, mesh_order=4): if not isinstance(resolution, (list, tuple)): resolution = (resolution,) * self.ambient_dim @@ -103,12 +125,24 @@ def get_mesh(self, resolution, mesh_order): order=mesh_order) +class BoxMeshBuilder1D(_BoxMeshBuilderBase): + ambient_dim = 1 + + +class BoxMeshBuilder2D(_BoxMeshBuilderBase): + ambient_dim = 2 + + +class BoxMeshBuilder3D(_BoxMeshBuilderBase): + ambient_dim = 2 + + class WarpedRectMeshBuilder(MeshBuilder): resolutions = [4, 6, 8] def __init__(self, dim): self.dim = dim - def get_mesh(self, resolution, mesh_order): + def get_mesh(self, resolution, mesh_order=4): return mgen.generate_warped_rect_mesh( dim=self.dim, order=4, nelements_side=6) diff --git a/test/test_dt_utils.py b/test/test_dt_utils.py index 9322f7d28..92d47dcdf 100644 --- a/test/test_dt_utils.py +++ b/test/test_dt_utils.py @@ -34,9 +34,10 @@ PytestPytatoPyOpenCLArrayContextFactory]) from grudge import DiscretizationCollection - import grudge.op as op +import mesh_data + import pytest import logging @@ -55,23 +56,22 @@ def test_geometric_factors_regular_refinement(actx_factory, name): # {{{ cases if name == "interval": - from mesh_data import BoxMeshBuilder - builder = BoxMeshBuilder(ambient_dim=1) + builder = mesh_data.BoxMeshBuilder1D() elif name == "box2d": - from mesh_data import BoxMeshBuilder - builder = BoxMeshBuilder(ambient_dim=2) + builder = mesh_data.BoxMeshBuilder2D() elif name == "box3d": - from mesh_data import BoxMeshBuilder - builder = BoxMeshBuilder(ambient_dim=3) + builder = mesh_data.BoxMeshBuilder3D() else: raise ValueError("unknown geometry name: %s" % name) # }}} + order = 4 + min_factors = [] for resolution in builder.resolutions: - mesh = builder.get_mesh(resolution, builder.mesh_order) - dcoll = DiscretizationCollection(actx, mesh, order=builder.order) + mesh = builder.get_mesh(resolution, order) + dcoll = DiscretizationCollection(actx, mesh, order=order) min_factors.append( actx.to_numpy( op.nodal_min(dcoll, "vol", actx.thaw(dt_geometric_factors(dcoll)))) @@ -84,8 +84,8 @@ def test_geometric_factors_regular_refinement(actx_factory, name): assert np.all(np.isclose(ratios, 2)) # Make sure it works with empty meshes - mesh = builder.get_mesh(0, builder.mesh_order) - dcoll = DiscretizationCollection(actx, mesh, order=builder.order) + mesh = builder.get_mesh(0) + dcoll = DiscretizationCollection(actx, mesh, order=order) factors = actx.thaw(dt_geometric_factors(dcoll)) # noqa: F841 @@ -98,14 +98,11 @@ def test_non_geometric_factors(actx_factory, name): # {{{ cases if name == "interval": - from mesh_data import BoxMeshBuilder - builder = BoxMeshBuilder(ambient_dim=1) + builder = mesh_data.BoxMeshBuilder1D() elif name == "box2d": - from mesh_data import BoxMeshBuilder - builder = BoxMeshBuilder(ambient_dim=2) + builder = mesh_data.BoxMeshBuilder2D() elif name == "box3d": - from mesh_data import BoxMeshBuilder - builder = BoxMeshBuilder(ambient_dim=3) + builder = mesh_data.BoxMeshBuilder3D() else: raise ValueError("unknown geometry name: %s" % name) diff --git a/test/test_grudge.py b/test/test_grudge.py index 3420ee95b..ff42eb63f 100644 --- a/test/test_grudge.py +++ b/test/test_grudge.py @@ -26,12 +26,15 @@ import numpy as np import numpy.linalg as la -from meshmode.discretization.poly_element import \ - InterpolatoryEdgeClusteredGroupFactory +from meshmode.discretization.poly_element import ( + InterpolatoryEdgeClusteredGroupFactory, + QuadratureGroupFactory) from meshmode.mesh import TensorProductElementGroup from grudge.array_context import PytestPyOpenCLArrayContextFactory from arraycontext import pytest_generate_tests_for_array_contexts +import mesh_data + pytest_generate_tests = pytest_generate_tests_for_array_contexts( [PytestPyOpenCLArrayContextFactory]) @@ -165,21 +168,19 @@ def test_mass_surface_area(actx_factory, name): # {{{ cases + order = 4 + if name == "2-1-ellipse": - from mesh_data import EllipseMeshBuilder - builder = EllipseMeshBuilder(radius=3.1, aspect_ratio=2.0) + builder = mesh_data.EllipseMeshBuilder(radius=3.1, aspect_ratio=2.0) surface_area = _ellipse_surface_area(builder.radius, builder.aspect_ratio) elif name == "spheroid": - from mesh_data import SpheroidMeshBuilder - builder = SpheroidMeshBuilder() + builder = mesh_data.SpheroidMeshBuilder() surface_area = _spheroid_surface_area(builder.radius, builder.aspect_ratio) elif name == "box2d": - from mesh_data import BoxMeshBuilder - builder = BoxMeshBuilder(ambient_dim=2) + builder = mesh_data.BoxMeshBuilder2D() surface_area = 1.0 elif name == "box3d": - from mesh_data import BoxMeshBuilder - builder = BoxMeshBuilder(ambient_dim=3) + builder = mesh_data.BoxMeshBuilder3D() surface_area = 1.0 else: raise ValueError("unknown geometry name: %s" % name) @@ -192,8 +193,8 @@ def test_mass_surface_area(actx_factory, name): eoc = EOCRecorder() for resolution in builder.resolutions: - mesh = builder.get_mesh(resolution, builder.mesh_order) - dcoll = DiscretizationCollection(actx, mesh, order=builder.order) + mesh = builder.get_mesh(resolution, order) + dcoll = DiscretizationCollection(actx, mesh, order=order) volume_discr = dcoll.discr_from_dd(dof_desc.DD_VOLUME) logger.info("ndofs: %d", volume_discr.ndofs) @@ -222,7 +223,7 @@ def test_mass_surface_area(actx_factory, name): logger.info("surface area error\n%s", str(eoc)) - assert eoc.max_error() < 3e-13 or eoc.order_estimate() > builder.order + assert eoc.max_error() < 3e-13 or eoc.order_estimate() > order # }}} @@ -234,13 +235,17 @@ def test_mass_surface_area(actx_factory, name): "spheroid", "warped_rect2", "warped_rect3", + "gh-339-1", + "gh-339-4", ]) def test_mass_operator_inverse(actx_factory, name): actx = actx_factory() # {{{ cases - import mesh_data + order = 4 + overintegrate = False + if name == "2-1-ellipse": # curve builder = mesh_data.EllipseMeshBuilder(radius=3.1, aspect_ratio=2.0) @@ -249,7 +254,19 @@ def test_mass_operator_inverse(actx_factory, name): builder = mesh_data.SpheroidMeshBuilder() elif name.startswith("warped_rect"): builder = mesh_data.WarpedRectMeshBuilder(dim=int(name[-1])) - + elif name == "gh-339-1": + builder = mesh_data.GmshMeshBuilder3D("gh-339.msh") + order = 1 + # NOTE: We're definitely not evaluating the bilinear forms accurately + # in that case, the mappings are very non-constant. + # It's kind of surprising that WADG manages to make a 15-digit inverse, + # but empirically it seems to. + elif name == "gh-339-1-overint": + builder = mesh_data.GmshMeshBuilder3D("gh-339.msh") + order = 1 + overintegrate = True + elif name == "gh-339-4": + builder = mesh_data.GmshMeshBuilder3D("gh-339.msh") else: raise ValueError("unknown geometry name: %s" % name) @@ -261,8 +278,14 @@ def test_mass_operator_inverse(actx_factory, name): eoc = EOCRecorder() for resolution in builder.resolutions: - mesh = builder.get_mesh(resolution, builder.mesh_order) - dcoll = DiscretizationCollection(actx, mesh, order=builder.order) + mesh = builder.get_mesh(resolution) + dcoll = make_discretization_collection( + actx, mesh, discr_tag_to_group_factory={ + dof_desc.DISCR_TAG_BASE: ( + InterpolatoryEdgeClusteredGroupFactory(order)), + dof_desc.DISCR_TAG_QUAD: ( + QuadratureGroupFactory(order)) + }) volume_discr = dcoll.discr_from_dd(dof_desc.DD_VOLUME) logger.info("ndofs: %d", volume_discr.ndofs) @@ -273,12 +296,21 @@ def test_mass_operator_inverse(actx_factory, name): def f(x): return actx.np.cos(4.0 * x[0]) - dd = dof_desc.DD_VOLUME x_volm = actx.thaw(volume_discr.nodes()) f_volm = f(x_volm) - f_inv = op.inverse_mass( - dcoll, op.mass(dcoll, dd, f_volm) - ) + if not overintegrate: + dd = dof_desc.DD_VOLUME + f_inv = op.inverse_mass( + dcoll, op.mass(dcoll, dd, f_volm) + ) + else: + dd_base_vol = dof_desc.as_dofdesc( + dof_desc.DTAG_VOLUME_ALL, dof_desc.DISCR_TAG_BASE) + dd_quad_vol = dof_desc.as_dofdesc( + dof_desc.DTAG_VOLUME_ALL, dof_desc.DISCR_TAG_QUAD) + f_inv = op.inverse_mass( + dcoll, op.mass(dcoll, dd_quad_vol, + op.project(dcoll, dd_base_vol, dd_quad_vol, f_volm))) inv_error = actx.to_numpy( op.norm(dcoll, f_volm - f_inv, 2) / op.norm(dcoll, f_volm, 2)) @@ -313,16 +345,15 @@ def test_face_normal_surface(actx_factory, mesh_name): # {{{ geometry if mesh_name == "2-1-ellipse": - from mesh_data import EllipseMeshBuilder - builder = EllipseMeshBuilder(radius=3.1, aspect_ratio=2.0) + builder = mesh_data.EllipseMeshBuilder(radius=3.1, aspect_ratio=2.0) elif mesh_name == "spheroid": - from mesh_data import SpheroidMeshBuilder - builder = SpheroidMeshBuilder() + builder = mesh_data.SpheroidMeshBuilder() else: raise ValueError("unknown mesh name: %s" % mesh_name) - mesh = builder.get_mesh(builder.resolutions[0], builder.mesh_order) - dcoll = DiscretizationCollection(actx, mesh, order=builder.order) + order = 4 + mesh = builder.get_mesh(builder.resolutions[0], order) + dcoll = DiscretizationCollection(actx, mesh, order=order) volume_discr = dcoll.discr_from_dd(dof_desc.DD_VOLUME) logger.info("ndofs: %d", volume_discr.ndofs) @@ -480,7 +511,6 @@ def test_gauss_theorem(actx_factory, case, visualize=False): raise ValueError(f"unknown case: {case}") from meshmode.mesh import BTAG_ALL - from meshmode.discretization.poly_element import QuadratureGroupFactory actx = actx_factory() @@ -572,20 +602,15 @@ def test_surface_divergence_theorem(actx_factory, mesh_name, visualize=False): # {{{ cases if mesh_name == "2-1-ellipse": - from mesh_data import EllipseMeshBuilder - builder = EllipseMeshBuilder(radius=3.1, aspect_ratio=2.0) + builder = mesh_data.EllipseMeshBuilder(radius=3.1, aspect_ratio=2.0) elif mesh_name == "2-1-spheroid": - from mesh_data import SpheroidMeshBuilder - builder = SpheroidMeshBuilder(radius=3.1, aspect_ratio=2.0) + builder = mesh_data.SpheroidMeshBuilder(radius=3.1, aspect_ratio=2.0) elif mesh_name == "circle": - from mesh_data import EllipseMeshBuilder - builder = EllipseMeshBuilder(radius=1.0, aspect_ratio=1.0) + builder = mesh_data.EllipseMeshBuilder(radius=1.0, aspect_ratio=1.0) elif mesh_name == "starfish": - from mesh_data import StarfishMeshBuilder - builder = StarfishMeshBuilder() + builder = mesh_data.StarfishMeshBuilder() elif mesh_name == "sphere": - from mesh_data import SphereMeshBuilder - builder = SphereMeshBuilder(radius=1.0, mesh_order=16) + builder = mesh_data.SphereMeshBuilder(radius=1.0) else: raise ValueError("unknown mesh name: %s" % mesh_name) @@ -618,13 +643,15 @@ def f(x): [0.0, np.sin(theta), np.cos(theta)], ]) + order = 4 + mesh_offset = np.array([0.33, -0.21, 0.0])[:ambient_dim] for i, resolution in enumerate(builder.resolutions): from meshmode.mesh.processing import affine_map from meshmode.discretization.connection import FACE_RESTR_ALL - mesh = builder.get_mesh(resolution, builder.mesh_order) + mesh = builder.get_mesh(resolution, order) mesh = affine_map(mesh, A=mesh_rotation, b=mesh_offset) from meshmode.discretization.poly_element import \ @@ -632,9 +659,9 @@ def f(x): qtag = dof_desc.DISCR_TAG_QUAD dcoll = DiscretizationCollection( - actx, mesh, order=builder.order, + actx, mesh, order=order, discr_tag_to_group_factory={ - qtag: QuadratureSimplexGroupFactory(2 * builder.order) + qtag: QuadratureSimplexGroupFactory(2 * order) } ) @@ -693,15 +720,15 @@ def f(x): # }}} - order = min(builder.order, builder.mesh_order) - 0.5 + exp_order = order - 0.5 logger.info("eoc_global:\n%s", str(eoc_global)) logger.info("eoc_local:\n%s", str(eoc_local)) assert eoc_global.max_error() < 1.0e-12 \ - or eoc_global.order_estimate() > order - 0.5 + or eoc_global.order_estimate() > exp_order - 0.5 assert eoc_local.max_error() < 1.0e-12 \ - or eoc_local.order_estimate() > order - 0.5 + or eoc_local.order_estimate() > exp_order - 0.5 # }}} diff --git a/test/test_reductions.py b/test/test_reductions.py index 3e4bf3964..ff858b47c 100644 --- a/test/test_reductions.py +++ b/test/test_reductions.py @@ -42,6 +42,7 @@ from grudge import DiscretizationCollection import grudge.op as op +import mesh_data from meshmode.dof_array import flatten @@ -64,11 +65,10 @@ def test_nodal_reductions(actx_factory, mesh_size, with_initial): actx = actx_factory() - from mesh_data import BoxMeshBuilder - builder = BoxMeshBuilder(ambient_dim=1) + builder = mesh_data.BoxMeshBuilder1D() - mesh = builder.get_mesh(mesh_size, builder.mesh_order) - dcoll = DiscretizationCollection(actx, mesh, order=builder.order) + mesh = builder.get_mesh(mesh_size) + dcoll = DiscretizationCollection(actx, mesh, order=4) x = actx.thaw(dcoll.nodes()) def f(x): @@ -130,12 +130,11 @@ def h(x): def test_elementwise_reductions(actx_factory): actx = actx_factory() - from mesh_data import BoxMeshBuilder - builder = BoxMeshBuilder(ambient_dim=1) + builder = mesh_data.BoxMeshBuilder1D() nelements = 4 - mesh = builder.get_mesh(nelements, builder.mesh_order) - dcoll = DiscretizationCollection(actx, mesh, order=builder.order) + mesh = builder.get_mesh(nelements) + dcoll = DiscretizationCollection(actx, mesh, order=4) x = actx.thaw(dcoll.nodes()) def f(x): @@ -192,11 +191,10 @@ def array_context(self): def test_nodal_reductions_with_container(actx_factory): actx = actx_factory() - from mesh_data import BoxMeshBuilder - builder = BoxMeshBuilder(ambient_dim=2) + builder = mesh_data.BoxMeshBuilder2D() - mesh = builder.get_mesh(4, builder.mesh_order) - dcoll = DiscretizationCollection(actx, mesh, order=builder.order) + mesh = builder.get_mesh(4) + dcoll = DiscretizationCollection(actx, mesh, order=4) x = actx.thaw(dcoll.nodes()) def f(x): @@ -239,12 +237,11 @@ def h(x): def test_elementwise_reductions_with_container(actx_factory): actx = actx_factory() - from mesh_data import BoxMeshBuilder - builder = BoxMeshBuilder(ambient_dim=2) + builder = mesh_data.BoxMeshBuilder2D() nelements = 4 - mesh = builder.get_mesh(nelements, builder.mesh_order) - dcoll = DiscretizationCollection(actx, mesh, order=builder.order) + mesh = builder.get_mesh(nelements) + dcoll = DiscretizationCollection(actx, mesh, order=4) x = actx.thaw(dcoll.nodes()) def f(x):