From b3eff8bddced86bc727cfd1dad471ffc5f344e3f Mon Sep 17 00:00:00 2001 From: Mike Campbell Date: Mon, 5 Feb 2024 16:09:50 -0600 Subject: [PATCH 1/5] Enable dt estimate for quads/hexes, extend tests for it Co-authored-by: Andreas Kloeckner --- grudge/dt_utils.py | 34 ++++++++++++++++++++++++---------- test/mesh_data.py | 2 ++ test/test_dt_utils.py | 36 ++++++++++++++++++++++++------------ 3 files changed, 50 insertions(+), 22 deletions(-) diff --git a/grudge/dt_utils.py b/grudge/dt_utils.py index 54a9343e8..0a9b223cf 100644 --- a/grudge/dt_utils.py +++ b/grudge/dt_utils.py @@ -229,20 +229,35 @@ def h_min_from_volume( def dt_geometric_factors( dcoll: DiscretizationCollection, dd: Optional[DOFDesc] = None) -> DOFArray: r"""Computes a geometric scaling factor for each cell following - [Hesthaven_2008]_, section 6.4, defined as the inradius (radius of an - inscribed circle/sphere). + [Hesthaven_2008]_, section 6.4, For simplicial elemenents, this factor is + defined as the inradius (radius of an inscribed circle/sphere). For + non-simplicial elements, a mean length measure is returned. - Specifically, the inradius for each element is computed using the following - formula from [Shewchuk_2002]_, Table 1, for simplicial cells - (triangles/tetrahedra): + Specifically, the inradius for each simplicial element is computed using the + following formula from [Shewchuk_2002]_, Table 1 (triangles, tetrahedra): .. math:: - r_D = \frac{d V}{\sum_{i=1}^{N_{faces}} F_i}, + r_D = \frac{d~V}{\sum_{i=1}^{N_{faces}} F_i}, where :math:`d` is the topological dimension, :math:`V` is the cell volume, and :math:`F_i` are the areas of each face of the cell. + For non-simplicial elements, we use the following formula for a mean + cell size measure: + + .. math:: + + r_D = \frac{2~d~V}{\sum_{i=1}^{N_{faces}} F_i}, + + where :math:`d` is the topological dimension, :math:`V` is the cell volume, + and :math:`F_i` are the areas of each face of the cell. Other valid choices + here include the shortest, longest, average of the cell diagonals, or edges. + The value returned by this routine (i.e. the cell volume divided by the + average cell face area) is bounded by the extrema of the cell edge lengths, + is straightforward to calculate regardless of element shape, and jibes well + with the foregoing calculation for simplicial elements. + :arg dd: a :class:`~grudge.dof_desc.DOFDesc`, or a value convertible to one. Defaults to the base volume discretization if not provided. :returns: a frozen :class:`~meshmode.dof_array.DOFArray` containing the @@ -256,11 +271,10 @@ def dt_geometric_factors( actx = dcoll._setup_actx volm_discr = dcoll.discr_from_dd(dd) + r_fac = dcoll.dim if any(not isinstance(grp, SimplexElementGroupBase) for grp in volm_discr.groups): - raise NotImplementedError( - "Geometric factors are only implemented for simplex element groups" - ) + r_fac = 2.0*r_fac if volm_discr.dim != volm_discr.ambient_dim: from warnings import warn @@ -342,7 +356,7 @@ def dt_geometric_factors( "e,ei->ei", 1/sae_i, actx.tag_axis(1, DiscretizationDOFAxisTag(), cv_i), - tagged=(FirstAxisIsElementsTag(),)) * dcoll.dim + tagged=(FirstAxisIsElementsTag(),)) * r_fac for cv_i, sae_i in zip(cell_vols, surface_areas))))) # }}} diff --git a/test/mesh_data.py b/test/mesh_data.py index 0ccc369a1..08799d09e 100644 --- a/test/mesh_data.py +++ b/test/mesh_data.py @@ -86,6 +86,7 @@ def get_mesh(self, resolution, mesh_order): class BoxMeshBuilder(MeshBuilder): ambient_dim = 2 + group_cls = None mesh_order = 1 resolutions = [4, 8, 16] @@ -100,6 +101,7 @@ def get_mesh(self, resolution, mesh_order): return mgen.generate_regular_rect_mesh( a=self.a, b=self.b, nelements_per_axis=resolution, + group_cls=self.group_cls, order=mesh_order) diff --git a/test/test_dt_utils.py b/test/test_dt_utils.py index 9322f7d28..84048db91 100644 --- a/test/test_dt_utils.py +++ b/test/test_dt_utils.py @@ -33,7 +33,7 @@ [PytestPyOpenCLArrayContextFactory, PytestPytatoPyOpenCLArrayContextFactory]) -from grudge import DiscretizationCollection +from grudge import make_discretization_collection import grudge.op as op @@ -47,22 +47,26 @@ @pytest.mark.parametrize("name", ["interval", "box2d", "box3d"]) -def test_geometric_factors_regular_refinement(actx_factory, name): +@pytest.mark.parametrize("tpe", [False, True]) +def test_geometric_factors_regular_refinement(actx_factory, name, tpe): from grudge.dt_utils import dt_geometric_factors actx = actx_factory() # {{{ cases + from meshmode.mesh import TensorProductElementGroup + group_cls = TensorProductElementGroup if tpe else None + if name == "interval": from mesh_data import BoxMeshBuilder - builder = BoxMeshBuilder(ambient_dim=1) + builder = BoxMeshBuilder(ambient_dim=1, group_cls=group_cls) elif name == "box2d": from mesh_data import BoxMeshBuilder - builder = BoxMeshBuilder(ambient_dim=2) + builder = BoxMeshBuilder(ambient_dim=2, group_cls=group_cls) elif name == "box3d": from mesh_data import BoxMeshBuilder - builder = BoxMeshBuilder(ambient_dim=3) + builder = BoxMeshBuilder(ambient_dim=3, group_cls=group_cls) else: raise ValueError("unknown geometry name: %s" % name) @@ -71,7 +75,7 @@ def test_geometric_factors_regular_refinement(actx_factory, name): min_factors = [] for resolution in builder.resolutions: mesh = builder.get_mesh(resolution, builder.mesh_order) - dcoll = DiscretizationCollection(actx, mesh, order=builder.order) + dcoll = make_discretization_collection(actx, mesh, order=builder.order) min_factors.append( actx.to_numpy( op.nodal_min(dcoll, "vol", actx.thaw(dt_geometric_factors(dcoll)))) @@ -85,7 +89,7 @@ def test_geometric_factors_regular_refinement(actx_factory, name): # Make sure it works with empty meshes mesh = builder.get_mesh(0, builder.mesh_order) - dcoll = DiscretizationCollection(actx, mesh, order=builder.order) + dcoll = make_discretization_collection(actx, mesh, order=builder.order) factors = actx.thaw(dt_geometric_factors(dcoll)) # noqa: F841 @@ -115,7 +119,7 @@ def test_non_geometric_factors(actx_factory, name): degrees = list(range(1, 8)) for degree in degrees: mesh = builder.get_mesh(1, degree) - dcoll = DiscretizationCollection(actx, mesh, order=degree) + dcoll = make_discretization_collection(actx, mesh, order=degree) factors.append(min(dt_non_geometric_factors(dcoll))) # Crude estimate, factors should behave like 1/N**2 @@ -134,7 +138,7 @@ def test_build_jacobian(actx_factory): mesh = mgen.generate_regular_rect_mesh(a=[0], b=[1], nelements_per_axis=(3,)) assert mesh.dim == 1 - dcoll = DiscretizationCollection(actx, mesh, order=1) + dcoll = make_discretization_collection(actx, mesh, order=1) def rhs(x): return 3*x**2 + 2*x + 5 @@ -151,19 +155,27 @@ def rhs(x): @pytest.mark.parametrize("dim", [1, 2]) @pytest.mark.parametrize("degree", [2, 4]) -def test_wave_dt_estimate(actx_factory, dim, degree, visualize=False): +@pytest.mark.parametrize("tpe", [False, True]) +def test_wave_dt_estimate(actx_factory, dim, degree, tpe, visualize=False): actx = actx_factory() + # {{{ cases + + from meshmode.mesh import TensorProductElementGroup + group_cls = TensorProductElementGroup if tpe else None + import meshmode.mesh.generation as mgen a = [0, 0, 0] b = [1, 1, 1] mesh = mgen.generate_regular_rect_mesh( a=a[:dim], b=b[:dim], - nelements_per_axis=(3,)*dim) + nelements_per_axis=(3,)*dim, + group_cls=group_cls) + assert mesh.dim == dim - dcoll = DiscretizationCollection(actx, mesh, order=degree) + dcoll = make_discretization_collection(actx, mesh, order=degree) from grudge.models.wave import WeakWaveOperator wave_op = WeakWaveOperator(dcoll, c=1) From 3406a9df16120328329366b09397d2fce293d941 Mon Sep 17 00:00:00 2001 From: Mike Campbell Date: Sat, 27 Jul 2024 01:06:33 -0500 Subject: [PATCH 2/5] Update dt_utils, wave, and tests for reasonable quad/hex DT estimates --- grudge/dt_utils.py | 133 +++++++++++++++++++++++++----------------- grudge/models/wave.py | 12 +++- test/mesh_data.py | 15 +++-- test/test_dt_utils.py | 39 ++++++------- 4 files changed, 115 insertions(+), 84 deletions(-) diff --git a/grudge/dt_utils.py b/grudge/dt_utils.py index a6b667ded..e5baf79f0 100644 --- a/grudge/dt_utils.py +++ b/grudge/dt_utils.py @@ -253,13 +253,13 @@ def dt_geometric_factors( .. math:: - r_D = \frac{2~d~V}{\sum_{i=1}^{N_{faces}} F_i}, + r_D = \frac{V}{2~\text{max}\left(F_i\right)}, where :math:`d` is the topological dimension, :math:`V` is the cell volume, and :math:`F_i` are the areas of each face of the cell. Other valid choices here include the shortest, longest, average of the cell diagonals, or edges. The value returned by this routine (i.e. the cell volume divided by the - average cell face area) is bounded by the extrema of the cell edge lengths, + maximum cell face area) is bounded by the extrema of the cell edge lengths, is straightforward to calculate regardless of element shape, and jibes well with the foregoing calculation for simplicial elements. @@ -276,10 +276,9 @@ def dt_geometric_factors( actx = dcoll._setup_actx volm_discr = dcoll.discr_from_dd(dd) - r_fac = dcoll.dim - if any(not isinstance(grp, SimplexElementGroupBase) - for grp in volm_discr.groups): - r_fac = 2.0*r_fac + tpe = any(not isinstance(grp, SimplexElementGroupBase) + for grp in volm_discr.groups) + r_fac = 0.5 if tpe else dcoll.dim if volm_discr.dim != volm_discr.ambient_dim: from warnings import warn @@ -308,61 +307,89 @@ def dt_geometric_factors( ) ) - if actx.supports_nonscalar_broadcasting: - # Compute total surface area of an element by summing over the - # individual face areas - surface_areas = DOFArray( - actx, - data=tuple( - actx.einsum( - "fej->e", - tag_axes(actx, { - 0: DiscretizationFaceAxisTag(), - 1: DiscretizationElementAxisTag(), - 2: DiscretizationDOFAxisTag() + if tpe: + if actx.supports_nonscalar_broadcasting: + surface_areas = DOFArray( + actx, + data=tuple( + actx.np.max( + tag_axes(actx, { + 0: DiscretizationFaceAxisTag(), + 1: DiscretizationElementAxisTag(), }, + face_ae_i.reshape( + vgrp.mesh_el_group.nfaces, + vgrp.nelements)), + axis=0) + for vgrp, face_ae_i in zip(volm_discr.groups, face_areas))) + else: + el_data_per_group = [] + for igrp, group in enumerate(volm_discr.mesh.groups): + nelements = group.nelements + nfaces = group.nfaces + el_face_data = face_areas[igrp].reshape( + nfaces, nelements, face_areas[igrp].shape[1]) + el_data_np = np.ascontiguousarray( + np.max(actx.to_numpy(el_face_data), axis=0)[:, 0:1]) + el_data = actx.from_numpy(el_data_np) + el_data = el_data.reshape(nelements) + el_data_per_group.append(el_data) + surface_areas = DOFArray(actx, tuple(el_data_per_group)) + else: + if actx.supports_nonscalar_broadcasting: + # Compute total surface area of an element by summing over the + # individual face areas + surface_areas = DOFArray( + actx, + data=tuple( + actx.einsum( + "fej->e", + tag_axes(actx, { + 0: DiscretizationFaceAxisTag(), + 1: DiscretizationElementAxisTag(), + 2: DiscretizationDOFAxisTag() + }, + face_ae_i.reshape( + vgrp.mesh_el_group.nfaces, + vgrp.nelements, + face_ae_i.shape[-1])), + tagged=(FirstAxisIsElementsTag(),)) + + for vgrp, face_ae_i in zip(volm_discr.groups, face_areas))) + else: + surface_areas = DOFArray( + actx, + data=tuple( + # NOTE: Whenever the array context can't perform nonscalar + # broadcasting, elementwise reductions + # (like `elementwise_integral`) repeat the *same* scalar value of + # the reduction at each degree of freedom. To get a single + # value for the face area (per face), + # we simply average over the nodes, which gives the desired result. + actx.einsum( + "fej->e", face_ae_i.reshape( vgrp.mesh_el_group.nfaces, vgrp.nelements, - face_ae_i.shape[-1])), - tagged=(FirstAxisIsElementsTag(),)) + face_ae_i.shape[-1] + ) / afgrp.nunit_dofs, + tagged=(FirstAxisIsElementsTag(),)) - for vgrp, face_ae_i in zip(volm_discr.groups, face_areas))) - else: - surface_areas = DOFArray( - actx, - data=tuple( - # NOTE: Whenever the array context can't perform nonscalar - # broadcasting, elementwise reductions - # (like `elementwise_integral`) repeat the *same* scalar value of - # the reduction at each degree of freedom. To get a single - # value for the face area (per face), - # we simply average over the nodes, which gives the desired result. - actx.einsum( - "fej->e", - face_ae_i.reshape( - vgrp.mesh_el_group.nfaces, - vgrp.nelements, - face_ae_i.shape[-1] - ) / afgrp.nunit_dofs, - tagged=(FirstAxisIsElementsTag(),)) - - for vgrp, afgrp, face_ae_i in zip(volm_discr.groups, - face_discr.groups, - face_areas) + for vgrp, afgrp, face_ae_i in zip(volm_discr.groups, + face_discr.groups, + face_areas) + ) ) - ) return actx.freeze( - actx.tag(NameHint(f"dt_geometric_{dd.as_identifier()}"), - DOFArray(actx, - data=tuple( - actx.einsum( - "e,ei->ei", - 1/sae_i, - actx.tag_axis(1, DiscretizationDOFAxisTag(), cv_i), - tagged=(FirstAxisIsElementsTag(),)) * r_fac - for cv_i, sae_i in zip(cell_vols, surface_areas))))) + actx.tag( + NameHint(f"dt_geometric_{dd.as_identifier()}"), + DOFArray(actx, data=tuple( + actx.einsum( + "e,ei->ei", 1/sae_i, + actx.tag_axis(1, DiscretizationDOFAxisTag(), cv_i), + tagged=(FirstAxisIsElementsTag(),)) * r_fac + for cv_i, sae_i in zip(cell_vols, surface_areas))))) # }}} diff --git a/grudge/models/wave.py b/grudge/models/wave.py index e1c645934..b5c6a8fa6 100644 --- a/grudge/models/wave.py +++ b/grudge/models/wave.py @@ -28,17 +28,18 @@ import numpy as np +from meshmode.discretization.poly_element import SimplexElementGroupBase from meshmode.mesh import BTAG_ALL, BTAG_NONE from pytools.obj_array import flat_obj_array import grudge.geometry as geo import grudge.op as op -from grudge.dof_desc import DISCR_TAG_BASE, as_dofdesc +from grudge.dof_desc import DISCR_TAG_BASE, as_dofdesc, DD_VOLUME_ALL from grudge.models import HyperbolicOperator - # {{{ constant-velocity + class WeakWaveOperator(HyperbolicOperator): r"""This operator discretizes the wave equation :math:`\partial_t^2 u = c^2 \Delta u`. @@ -187,7 +188,12 @@ def max_characteristic_velocity(self, actx, t=None, fields=None): def estimate_rk4_timestep(self, actx, dcoll, **kwargs): # FIXME: Sketchy, empirically determined fudge factor - return 0.38 * super().estimate_rk4_timestep(actx, dcoll, **kwargs) + volm_discr = dcoll.discr_from_dd(DD_VOLUME_ALL) + tpe = any(not isinstance(grp, SimplexElementGroupBase) + for grp in volm_discr.groups) + fudge_factor = 0.23 if tpe else 0.38 + return fudge_factor * super().estimate_rk4_timestep( + actx, dcoll, **kwargs) # }}} diff --git a/test/mesh_data.py b/test/mesh_data.py index b331a7449..cb3435242 100644 --- a/test/mesh_data.py +++ b/test/mesh_data.py @@ -4,13 +4,14 @@ import numpy as np import meshmode.mesh.generation as mgen -from meshmode.mesh import Mesh +from meshmode.mesh import Mesh, MeshElementGroup from meshmode.mesh.io import read_gmsh class MeshBuilder(ABC): resolutions: ClassVar[Sequence[Hashable]] ambient_dim: ClassVar[int] + group_cls: ClassVar[MeshElementGroup] @abstractmethod def get_mesh( @@ -111,8 +112,10 @@ def get_mesh(self, resolution, mesh_order=4): class BoxMeshBuilder(MeshBuilder): - ambient_dim = 2 - group_cls = None + + def __init__(self, ambient_dim=2, group_cls=None): + self.ambient_dim = ambient_dim + self.group_cls = group_cls mesh_order = 1 @@ -130,15 +133,15 @@ def get_mesh(self, resolution, mesh_order=4): order=mesh_order) -class BoxMeshBuilder1D(_BoxMeshBuilderBase): +class BoxMeshBuilder1D(BoxMeshBuilder): ambient_dim = 1 -class BoxMeshBuilder2D(_BoxMeshBuilderBase): +class BoxMeshBuilder2D(BoxMeshBuilder): ambient_dim = 2 -class BoxMeshBuilder3D(_BoxMeshBuilderBase): +class BoxMeshBuilder3D(BoxMeshBuilder): ambient_dim = 2 diff --git a/test/test_dt_utils.py b/test/test_dt_utils.py index 3df99e509..63d2b6e3c 100644 --- a/test/test_dt_utils.py +++ b/test/test_dt_utils.py @@ -30,26 +30,18 @@ PytestPyOpenCLArrayContextFactory, PytestPytatoPyOpenCLArrayContextFactory, ) - - -pytest_generate_tests = pytest_generate_tests_for_array_contexts( - [PytestPyOpenCLArrayContextFactory, - PytestPytatoPyOpenCLArrayContextFactory]) - -from grudge import make_discretization_collection - +from grudge.discretization import make_discretization_collection import grudge.op as op import mesh_data import pytest -import grudge.op as op -from grudge.discretization import make_discretization_collection - - -logger = logging.getLogger(__name__) from meshmode import _acf # noqa: F401 +pytest_generate_tests = pytest_generate_tests_for_array_contexts( + [PytestPyOpenCLArrayContextFactory, + PytestPytatoPyOpenCLArrayContextFactory]) + @pytest.mark.parametrize("name", ["interval", "box2d", "box3d"]) @pytest.mark.parametrize("tpe", [False, True]) @@ -78,14 +70,15 @@ def test_geometric_factors_regular_refinement(actx_factory, name, tpe): # }}} order = 4 - min_factors = [] - for resolution in builder.resolutions: - mesh = builder.get_mesh(resolution, builder.mesh_order) - dcoll = make_discretization_collection(actx, mesh, order=builder.order) + resolutions = [2, 4, 8] + for resolution in resolutions: + mesh = builder.get_mesh(resolution, order) + dcoll = make_discretization_collection(actx, mesh, order=order) min_factors.append( actx.to_numpy( - op.nodal_min(dcoll, "vol", actx.thaw(dt_geometric_factors(dcoll)))) + op.nodal_min(dcoll, "vol", + actx.thaw(dt_geometric_factors(dcoll)))) ) # Resolution is doubled each refinement, so the ratio of consecutive @@ -95,9 +88,10 @@ def test_geometric_factors_regular_refinement(actx_factory, name, tpe): assert np.all(np.isclose(ratios, 2)) # Make sure it works with empty meshes - mesh = builder.get_mesh(0, builder.mesh_order) - dcoll = make_discretization_collection(actx, mesh, order=builder.order) - factors = actx.thaw(dt_geometric_factors(dcoll)) # noqa: F841 + if not tpe: + mesh = builder.get_mesh(0, order) + dcoll = make_discretization_collection(actx, mesh, order=order) + factors = actx.thaw(dt_geometric_factors(dcoll)) # noqa: F841 @pytest.mark.parametrize("name", ["interval", "box2d", "box3d"]) @@ -164,10 +158,11 @@ def test_wave_dt_estimate(actx_factory, dim, degree, tpe, visualize=False): actx = actx_factory() # {{{ cases + if dim == 1 and tpe: + pytest.skip("Skipping 1D test for tensor product elements.") from meshmode.mesh import TensorProductElementGroup group_cls = TensorProductElementGroup if tpe else None - import meshmode.mesh.generation as mgen a = [0, 0, 0] From 24e17aeb5c55a2866442dfd4a4086b49d63c5ea6 Mon Sep 17 00:00:00 2001 From: Mike Campbell Date: Sat, 27 Jul 2024 07:47:59 -0500 Subject: [PATCH 3/5] Scooch closer to main --- grudge/models/wave.py | 4 ++-- test/mesh_data.py | 1 + test/test_dt_utils.py | 21 +++++++++++++-------- 3 files changed, 16 insertions(+), 10 deletions(-) diff --git a/grudge/models/wave.py b/grudge/models/wave.py index b5c6a8fa6..b22ff1594 100644 --- a/grudge/models/wave.py +++ b/grudge/models/wave.py @@ -34,11 +34,11 @@ import grudge.geometry as geo import grudge.op as op -from grudge.dof_desc import DISCR_TAG_BASE, as_dofdesc, DD_VOLUME_ALL +from grudge.dof_desc import DD_VOLUME_ALL, DISCR_TAG_BASE, as_dofdesc from grudge.models import HyperbolicOperator -# {{{ constant-velocity +# {{{ constant-velocity class WeakWaveOperator(HyperbolicOperator): r"""This operator discretizes the wave equation diff --git a/test/mesh_data.py b/test/mesh_data.py index cb3435242..dff5cfa14 100644 --- a/test/mesh_data.py +++ b/test/mesh_data.py @@ -112,6 +112,7 @@ def get_mesh(self, resolution, mesh_order=4): class BoxMeshBuilder(MeshBuilder): + resolutions: ClassVar[Sequence[Hashable]] = [4, 8, 16] def __init__(self, ambient_dim=2, group_cls=None): self.ambient_dim = ambient_dim diff --git a/test/test_dt_utils.py b/test/test_dt_utils.py index 63d2b6e3c..aa6ebc099 100644 --- a/test/test_dt_utils.py +++ b/test/test_dt_utils.py @@ -30,18 +30,24 @@ PytestPyOpenCLArrayContextFactory, PytestPytatoPyOpenCLArrayContextFactory, ) -from grudge.discretization import make_discretization_collection -import grudge.op as op - -import mesh_data -import pytest -from meshmode import _acf # noqa: F401 pytest_generate_tests = pytest_generate_tests_for_array_contexts( [PytestPyOpenCLArrayContextFactory, PytestPytatoPyOpenCLArrayContextFactory]) +import logging + +import mesh_data +import pytest + +import grudge.op as op +from grudge.discretization import make_discretization_collection + + +logger = logging.getLogger(__name__) +from meshmode import _acf # noqa: F401 + @pytest.mark.parametrize("name", ["interval", "box2d", "box3d"]) @pytest.mark.parametrize("tpe", [False, True]) @@ -71,8 +77,7 @@ def test_geometric_factors_regular_refinement(actx_factory, name, tpe): order = 4 min_factors = [] - resolutions = [2, 4, 8] - for resolution in resolutions: + for resolution in builder.resolutions: mesh = builder.get_mesh(resolution, order) dcoll = make_discretization_collection(actx, mesh, order=order) min_factors.append( From 7b8f37b4547873d87baf2444405a583c5100fe97 Mon Sep 17 00:00:00 2001 From: Mike Campbell Date: Sun, 28 Jul 2024 14:35:05 -0500 Subject: [PATCH 4/5] Include warpy mesh in dt test. --- grudge/models/wave.py | 2 +- test/test_dt_utils.py | 49 +++++++++++++++++++++++++++---------------- 2 files changed, 32 insertions(+), 19 deletions(-) diff --git a/grudge/models/wave.py b/grudge/models/wave.py index b22ff1594..9f07797d9 100644 --- a/grudge/models/wave.py +++ b/grudge/models/wave.py @@ -191,7 +191,7 @@ def estimate_rk4_timestep(self, actx, dcoll, **kwargs): volm_discr = dcoll.discr_from_dd(DD_VOLUME_ALL) tpe = any(not isinstance(grp, SimplexElementGroupBase) for grp in volm_discr.groups) - fudge_factor = 0.23 if tpe else 0.38 + fudge_factor = 0.233 if tpe else 0.38 return fudge_factor * super().estimate_rk4_timestep( actx, dcoll, **kwargs) diff --git a/test/test_dt_utils.py b/test/test_dt_utils.py index aa6ebc099..04375a37f 100644 --- a/test/test_dt_utils.py +++ b/test/test_dt_utils.py @@ -47,6 +47,7 @@ logger = logging.getLogger(__name__) from meshmode import _acf # noqa: F401 +from meshmode.mesh.io import read_gmsh @pytest.mark.parametrize("name", ["interval", "box2d", "box3d"]) @@ -156,35 +157,46 @@ def rhs(x): assert np.allclose(np.diag(mat), 3*2*2 + 2) -@pytest.mark.parametrize("dim", [1, 2]) +# ("simple_tets.msh", 3, False, False), +# ("simple_hexs.msh", 3, True, False)]) @pytest.mark.parametrize("degree", [2, 4]) -@pytest.mark.parametrize("tpe", [False, True]) -def test_wave_dt_estimate(actx_factory, dim, degree, tpe, visualize=False): +@pytest.mark.parametrize(("meshfile", "dim", "tpe", "warp"), + [(None, 1, False, False), + (None, 2, False, False), + (None, 2, False, True), + (None, 2, True, False), + (None, 2, True, True),]) +def test_wave_dt_estimate(actx_factory, degree, meshfile, dim, tpe, warp, + visualize=False): actx = actx_factory() # {{{ cases - if dim == 1 and tpe: - pytest.skip("Skipping 1D test for tensor product elements.") - from meshmode.mesh import TensorProductElementGroup group_cls = TensorProductElementGroup if tpe else None import meshmode.mesh.generation as mgen - a = [0, 0, 0] - b = [1, 1, 1] - mesh = mgen.generate_regular_rect_mesh( - a=a[:dim], b=b[:dim], - nelements_per_axis=(3,)*dim, - group_cls=group_cls) - - assert mesh.dim == dim - + if meshfile is None: + if warp: + mesh = mgen.generate_warped_rect_mesh( + dim=dim, order=2, nelements_side=3, group_cls=group_cls) + else: + a = [0, 0, 0] + b = [1, 1, 1] + mesh = mgen.generate_regular_rect_mesh( + a=a[:dim], b=b[:dim], + nelements_per_axis=(3,)*dim, + group_cls=group_cls) + assert mesh.dim == dim + else: + mesh_construction_kwargs = {"force_positive_orientation": True} + mesh = read_gmsh(meshfile, + mesh_construction_kwargs=mesh_construction_kwargs) dcoll = make_discretization_collection(actx, mesh, order=degree) from grudge.models.wave import WeakWaveOperator wave_op = WeakWaveOperator(dcoll, c=1) rhs = actx.compile( - lambda w: wave_op.operator(t=0, w=w)) + lambda w: wave_op.operator(t=0, w=w)) from pytools.obj_array import make_obj_array fields = make_obj_array([dcoll.zeros(actx) for i in range(dim+1)]) @@ -212,7 +224,8 @@ def test_wave_dt_estimate(actx_factory, dim, degree, tpe, visualize=False): import matplotlib.pyplot as plt plt.contour(re, im, sf_grid, [0.25, 0.5, 0.75, 0.9, 1, 1.1]) plt.colorbar() - plt.plot(dt_est * eigvals.real, dt_est * eigvals.imag, "x") + plt.plot(dt_est * eigvals.real, dt_est * eigvals.imag, + "x") plt.grid() plt.show() @@ -230,7 +243,7 @@ def test_wave_dt_estimate(actx_factory, dim, degree, tpe, visualize=False): print(f"Stable timestep is {max(stable_dt_factors):.2f}x the estimate") else: print("Stable timestep estimate appears to be sharp") - assert not stable_dt_factors or max(stable_dt_factors) < 1.5, stable_dt_factors + assert not stable_dt_factors or max(stable_dt_factors) < 1.54, stable_dt_factors # You can test individual routines by typing From 3c4ebeae9feba41cd91ee5b276a9abfd2cd24e3b Mon Sep 17 00:00:00 2001 From: Mike Campbell Date: Sun, 28 Jul 2024 14:40:10 -0500 Subject: [PATCH 5/5] Add some simple gmsh test meshes. --- test/simple_hexs.msh | 356 +++++++++++++++++++++++++++++++++++++++++++ test/simple_tets.msh | 91 +++++++++++ 2 files changed, 447 insertions(+) create mode 100644 test/simple_hexs.msh create mode 100644 test/simple_tets.msh diff --git a/test/simple_hexs.msh b/test/simple_hexs.msh new file mode 100644 index 000000000..1b271835b --- /dev/null +++ b/test/simple_hexs.msh @@ -0,0 +1,356 @@ +$MeshFormat +2.2 0 8 +$EndMeshFormat +$Nodes +147 +1 0 0 1 +2 0 0 0 +3 0 1 1 +4 0 1 0 +5 1 0 1 +6 1 0 0 +7 1 1 1 +8 1 1 0 +9 0 0 0.4999999999999999 +10 0 0.4999999999999999 1 +11 0 1 0.4999999999999999 +12 0 0.4999999999999999 0 +13 1 0 0.4999999999999999 +14 1 0.4999999999999999 1 +15 1 1 0.4999999999999999 +16 1 0.4999999999999999 0 +17 0.4999999999999999 0 0 +18 0.4999999999999999 0 1 +19 0.4999999999999999 1 0 +20 0.4999999999999999 1 1 +21 0 0.5 0.5 +22 0 0.25 0.75 +23 0 0.25 0.25 +24 0 0.75 0.75 +25 0 0.75 0.25 +26 0 0.1666666666666667 0.5 +27 0 0.5 0.8333333333333334 +28 0 0.5 0.1666666666666667 +29 0 0.8333333333333334 0.5 +30 1 0.5 0.5 +31 1 0.25 0.25 +32 1 0.25 0.75 +33 1 0.75 0.75 +34 1 0.75 0.25 +35 1 0.1666666666666667 0.5 +36 1 0.5 0.8333333333333334 +37 1 0.5 0.1666666666666667 +38 1 0.8333333333333334 0.5 +39 0.5 0 0.5 +40 0.25 0 0.25 +41 0.25 0 0.75 +42 0.75 0 0.75 +43 0.75 0 0.25 +44 0.1666666666666667 0 0.5 +45 0.5 0 0.8333333333333334 +46 0.5 0 0.1666666666666667 +47 0.8333333333333334 0 0.5 +48 0.5 1 0.5 +49 0.25 1 0.75 +50 0.25 1 0.25 +51 0.75 1 0.75 +52 0.75 1 0.25 +53 0.1666666666666667 1 0.5 +54 0.5 1 0.8333333333333334 +55 0.5 1 0.1666666666666667 +56 0.8333333333333334 1 0.5 +57 0.5 0.5 0 +58 0.25 0.75 0 +59 0.25 0.25 0 +60 0.75 0.25 0 +61 0.75 0.75 0 +62 0.1666666666666667 0.5 0 +63 0.5 0.1666666666666667 0 +64 0.5 0.8333333333333334 0 +65 0.8333333333333334 0.5 0 +66 0.5 0.5 1 +67 0.25 0.25 1 +68 0.25 0.75 1 +69 0.75 0.25 1 +70 0.75 0.75 1 +71 0.1666666666666667 0.5 1 +72 0.5 0.1666666666666667 1 +73 0.5 0.8333333333333334 1 +74 0.8333333333333334 0.5 1 +75 0.75 0.5 0.25 +76 0.5 0.5 0.5 +77 0.25 0.5 0.25 +78 0.5 0.25 0.25 +79 0.25 0.25 0.5 +80 0.75 0.25 0.5 +81 0.25 0.5 0.75 +82 0.25 0.75 0.5 +83 0.5 0.75 0.75 +84 0.75 0.5 0.75 +85 0.75 0.75 0.5 +86 0.5 0.75 0.25 +87 0.5 0.25 0.75 +88 0.5 0.5 0.3333333333333333 +89 0.6666666666666666 0.3333333333333333 0.3333333333333333 +90 0.3333333333333333 0.3333333333333333 0.3333333333333333 +91 0.5 0.3333333333333333 0.5 +92 0.5 0.375 0.375 +93 0.3333333333333333 0.6666666666666666 0.6666666666666666 +94 0.5 0.5 0.6666666666666666 +95 0.6666666666666666 0.6666666666666666 0.6666666666666666 +96 0.5 0.6666666666666666 0.5 +97 0.5 0.625 0.625 +98 0.6666666666666666 0.6666666666666666 0.3333333333333333 +99 0.3333333333333333 0.6666666666666666 0.3333333333333333 +100 0.5 0.625 0.375 +101 0.3333333333333333 0.3333333333333333 0.6666666666666666 +102 0.6666666666666666 0.3333333333333333 0.6666666666666666 +103 0.5 0.375 0.625 +104 0.8333333333333334 0.8333333333333334 0.6666666666666666 +105 0.8333333333333334 0.8333333333333334 0.3333333333333333 +106 0.875 0.875 0.5 +107 0.1666666666666667 0.8333333333333334 0.3333333333333333 +108 0.1666666666666667 0.8333333333333334 0.6666666666666666 +109 0.125 0.875 0.5 +110 0.6666666666666666 0.1666666666666667 0.8333333333333334 +111 0.3333333333333333 0.1666666666666667 0.8333333333333334 +112 0.5 0.125 0.875 +113 0.1666666666666667 0.6666666666666666 0.8333333333333334 +114 0.1666666666666667 0.3333333333333333 0.8333333333333334 +115 0.125 0.5 0.875 +116 0.3333333333333333 0.8333333333333334 0.8333333333333334 +117 0.6666666666666666 0.8333333333333334 0.8333333333333334 +118 0.5 0.875 0.875 +119 0.8333333333333334 0.6666666666666666 0.1666666666666667 +120 0.8333333333333334 0.3333333333333333 0.1666666666666667 +121 0.875 0.5 0.125 +122 0.8333333333333334 0.6666666666666666 0.8333333333333334 +123 0.8333333333333334 0.3333333333333333 0.8333333333333334 +124 0.875 0.5 0.875 +125 0.1666666666666667 0.3333333333333333 0.1666666666666667 +126 0.1666666666666667 0.6666666666666666 0.1666666666666667 +127 0.125 0.5 0.125 +128 0.1666666666666667 0.1666666666666667 0.6666666666666666 +129 0.1666666666666667 0.1666666666666667 0.3333333333333333 +130 0.125 0.125 0.5 +131 0.3333333333333333 0.8333333333333334 0.1666666666666667 +132 0.6666666666666666 0.8333333333333334 0.1666666666666667 +133 0.5 0.875 0.125 +134 0.6666666666666666 0.1666666666666667 0.1666666666666667 +135 0.3333333333333333 0.1666666666666667 0.1666666666666667 +136 0.5 0.125 0.125 +137 0.8333333333333334 0.1666666666666667 0.3333333333333333 +138 0.8333333333333334 0.1666666666666667 0.6666666666666666 +139 0.875 0.125 0.5 +140 0.25 0.25 0.75 +141 0.25 0.75 0.75 +142 0.75 0.75 0.75 +143 0.25 0.75 0.25 +144 0.75 0.75 0.25 +145 0.75 0.25 0.25 +146 0.75 0.25 0.75 +147 0.25 0.25 0.25 +$EndNodes +$Elements +200 +1 15 2 0 1 1 +2 15 2 0 2 2 +3 15 2 0 3 3 +4 15 2 0 4 4 +5 15 2 0 5 5 +6 15 2 0 6 6 +7 15 2 0 7 7 +8 15 2 0 8 8 +9 1 2 0 1 2 9 +10 1 2 0 1 9 1 +11 1 2 0 2 1 10 +12 1 2 0 2 10 3 +13 1 2 0 3 4 11 +14 1 2 0 3 11 3 +15 1 2 0 4 2 12 +16 1 2 0 4 12 4 +17 1 2 0 5 6 13 +18 1 2 0 5 13 5 +19 1 2 0 6 5 14 +20 1 2 0 6 14 7 +21 1 2 0 7 8 15 +22 1 2 0 7 15 7 +23 1 2 0 8 6 16 +24 1 2 0 8 16 8 +25 1 2 0 9 2 17 +26 1 2 0 9 17 6 +27 1 2 0 10 1 18 +28 1 2 0 10 18 5 +29 1 2 0 11 4 19 +30 1 2 0 11 19 8 +31 1 2 0 12 3 20 +32 1 2 0 12 20 7 +33 3 2 0 1 2 9 26 23 +34 3 2 0 1 9 1 22 26 +35 3 2 0 1 23 26 22 21 +36 3 2 0 1 1 10 27 22 +37 3 2 0 1 10 3 24 27 +38 3 2 0 1 22 27 24 21 +39 3 2 0 1 4 12 28 25 +40 3 2 0 1 12 2 23 28 +41 3 2 0 1 25 28 23 21 +42 3 2 0 1 3 11 29 24 +43 3 2 0 1 11 4 25 29 +44 3 2 0 1 24 29 25 21 +45 3 2 0 2 6 31 35 13 +46 3 2 0 2 31 30 32 35 +47 3 2 0 2 13 35 32 5 +48 3 2 0 2 5 32 36 14 +49 3 2 0 2 32 30 33 36 +50 3 2 0 2 14 36 33 7 +51 3 2 0 2 8 34 37 16 +52 3 2 0 2 34 30 31 37 +53 3 2 0 2 16 37 31 6 +54 3 2 0 2 7 33 38 15 +55 3 2 0 2 33 30 34 38 +56 3 2 0 2 15 38 34 8 +57 3 2 0 3 1 9 44 41 +58 3 2 0 3 9 2 40 44 +59 3 2 0 3 41 44 40 39 +60 3 2 0 3 5 18 45 42 +61 3 2 0 3 18 1 41 45 +62 3 2 0 3 42 45 41 39 +63 3 2 0 3 2 17 46 40 +64 3 2 0 3 17 6 43 46 +65 3 2 0 3 40 46 43 39 +66 3 2 0 3 6 13 47 43 +67 3 2 0 3 13 5 42 47 +68 3 2 0 3 43 47 42 39 +69 3 2 0 4 3 49 53 11 +70 3 2 0 4 49 48 50 53 +71 3 2 0 4 11 53 50 4 +72 3 2 0 4 7 51 54 20 +73 3 2 0 4 51 48 49 54 +74 3 2 0 4 20 54 49 3 +75 3 2 0 4 4 50 55 19 +76 3 2 0 4 50 48 52 55 +77 3 2 0 4 19 55 52 8 +78 3 2 0 4 8 52 56 15 +79 3 2 0 4 52 48 51 56 +80 3 2 0 4 15 56 51 7 +81 3 2 0 5 2 12 62 59 +82 3 2 0 5 12 4 58 62 +83 3 2 0 5 59 62 58 57 +84 3 2 0 5 6 17 63 60 +85 3 2 0 5 17 2 59 63 +86 3 2 0 5 60 63 59 57 +87 3 2 0 5 4 19 64 58 +88 3 2 0 5 19 8 61 64 +89 3 2 0 5 58 64 61 57 +90 3 2 0 5 8 16 65 61 +91 3 2 0 5 16 6 60 65 +92 3 2 0 5 61 65 60 57 +93 3 2 0 6 1 67 71 10 +94 3 2 0 6 67 66 68 71 +95 3 2 0 6 10 71 68 3 +96 3 2 0 6 5 69 72 18 +97 3 2 0 6 69 66 67 72 +98 3 2 0 6 18 72 67 1 +99 3 2 0 6 3 68 73 20 +100 3 2 0 6 68 66 70 73 +101 3 2 0 6 20 73 70 7 +102 3 2 0 6 7 70 74 14 +103 3 2 0 6 70 66 69 74 +104 3 2 0 6 14 74 69 5 +105 5 2 0 1 57 75 88 77 78 89 92 90 +106 5 2 0 1 75 30 76 88 89 80 91 92 +107 5 2 0 1 77 88 76 21 90 92 91 79 +108 5 2 0 1 39 80 89 78 79 91 92 90 +109 5 2 0 1 66 81 93 83 84 94 97 95 +110 5 2 0 1 81 21 82 93 94 76 96 97 +111 5 2 0 1 83 93 82 48 95 97 96 85 +112 5 2 0 1 30 76 94 84 85 96 97 95 +113 5 2 0 1 30 76 96 85 75 88 100 98 +114 5 2 0 1 76 21 82 96 88 77 99 100 +115 5 2 0 1 85 96 82 48 98 100 99 86 +116 5 2 0 1 57 77 88 75 86 99 100 98 +117 5 2 0 1 66 81 94 84 87 101 103 102 +118 5 2 0 1 81 21 76 94 101 79 91 103 +119 5 2 0 1 84 94 76 30 102 103 91 80 +120 5 2 0 1 39 79 101 87 80 91 103 102 +121 5 2 0 1 48 51 104 85 52 56 106 105 +122 5 2 0 1 51 7 33 104 56 15 38 106 +123 5 2 0 1 85 104 33 30 105 106 38 34 +124 5 2 0 1 8 15 56 52 34 38 106 105 +125 5 2 0 1 21 82 107 25 24 108 109 29 +126 5 2 0 1 82 48 50 107 108 49 53 109 +127 5 2 0 1 25 107 50 4 29 109 53 11 +128 5 2 0 1 3 49 108 24 11 53 109 29 +129 5 2 0 1 5 18 72 69 42 45 112 110 +130 5 2 0 1 18 1 67 72 45 41 111 112 +131 5 2 0 1 69 72 67 66 110 112 111 87 +132 5 2 0 1 39 41 45 42 87 111 112 110 +133 5 2 0 1 66 68 113 81 67 71 115 114 +134 5 2 0 1 68 3 24 113 71 10 27 115 +135 5 2 0 1 81 113 24 21 114 115 27 22 +136 5 2 0 1 1 10 71 67 22 27 115 114 +137 5 2 0 1 48 49 116 83 51 54 118 117 +138 5 2 0 1 49 3 68 116 54 20 73 118 +139 5 2 0 1 83 116 68 66 117 118 73 70 +140 5 2 0 1 7 20 54 51 70 73 118 117 +141 5 2 0 1 57 61 119 75 60 65 121 120 +142 5 2 0 1 61 8 34 119 65 16 37 121 +143 5 2 0 1 75 119 34 30 120 121 37 31 +144 5 2 0 1 6 16 65 60 31 37 121 120 +145 5 2 0 1 7 14 74 70 33 36 124 122 +146 5 2 0 1 14 5 69 74 36 32 123 124 +147 5 2 0 1 70 74 69 66 122 124 123 84 +148 5 2 0 1 30 32 36 33 84 123 124 122 +149 5 2 0 1 21 23 28 25 77 125 127 126 +150 5 2 0 1 23 2 12 28 125 59 62 127 +151 5 2 0 1 25 28 12 4 126 127 62 58 +152 5 2 0 1 57 59 125 77 58 62 127 126 +153 5 2 0 1 1 9 26 22 41 44 130 128 +154 5 2 0 1 9 2 23 26 44 40 129 130 +155 5 2 0 1 22 26 23 21 128 130 129 79 +156 5 2 0 1 39 40 44 41 79 129 130 128 +157 5 2 0 1 48 50 55 52 86 131 133 132 +158 5 2 0 1 50 4 19 55 131 58 64 133 +159 5 2 0 1 52 55 19 8 132 133 64 61 +160 5 2 0 1 57 58 131 86 61 64 133 132 +161 5 2 0 1 6 60 63 17 43 134 136 46 +162 5 2 0 1 60 57 59 63 134 78 135 136 +163 5 2 0 1 17 63 59 2 46 136 135 40 +164 5 2 0 1 39 78 134 43 40 135 136 46 +165 5 2 0 1 6 13 35 31 43 47 139 137 +166 5 2 0 1 13 5 32 35 47 42 138 139 +167 5 2 0 1 31 35 32 30 137 139 138 80 +168 5 2 0 1 39 42 47 43 80 138 139 137 +169 5 2 0 1 21 81 114 22 79 101 140 128 +170 5 2 0 1 81 66 67 114 101 87 111 140 +171 5 2 0 1 22 114 67 1 128 140 111 41 +172 5 2 0 1 39 87 101 79 41 111 140 128 +173 5 2 0 1 3 24 108 49 68 113 141 116 +174 5 2 0 1 24 21 82 108 113 81 93 141 +175 5 2 0 1 49 108 82 48 116 141 93 83 +176 5 2 0 1 66 81 113 68 83 93 141 116 +177 5 2 0 1 66 83 117 70 84 95 142 122 +178 5 2 0 1 83 48 51 117 95 85 104 142 +179 5 2 0 1 70 117 51 7 122 142 104 33 +180 5 2 0 1 30 85 95 84 33 104 142 122 +181 5 2 0 1 4 50 107 25 58 131 143 126 +182 5 2 0 1 50 48 82 107 131 86 99 143 +183 5 2 0 1 25 107 82 21 126 143 99 77 +184 5 2 0 1 57 86 131 58 77 99 143 126 +185 5 2 0 1 8 34 105 52 61 119 144 132 +186 5 2 0 1 34 30 85 105 119 75 98 144 +187 5 2 0 1 52 105 85 48 132 144 98 86 +188 5 2 0 1 57 75 119 61 86 98 144 132 +189 5 2 0 1 30 75 120 31 80 89 145 137 +190 5 2 0 1 75 57 60 120 89 78 134 145 +191 5 2 0 1 31 120 60 6 137 145 134 43 +192 5 2 0 1 39 78 89 80 43 134 145 137 +193 5 2 0 1 30 32 123 84 80 138 146 102 +194 5 2 0 1 32 5 69 123 138 42 110 146 +195 5 2 0 1 84 123 69 66 102 146 110 87 +196 5 2 0 1 39 42 138 80 87 110 146 102 +197 5 2 0 1 21 23 125 77 79 129 147 90 +198 5 2 0 1 23 2 59 125 129 40 135 147 +199 5 2 0 1 77 125 59 57 90 147 135 78 +200 5 2 0 1 39 40 129 79 78 135 147 90 +$EndElements diff --git a/test/simple_tets.msh b/test/simple_tets.msh new file mode 100644 index 000000000..eb1cd7dca --- /dev/null +++ b/test/simple_tets.msh @@ -0,0 +1,91 @@ +$MeshFormat +2.2 0 8 +$EndMeshFormat +$Nodes +14 +1 0 0 1 +2 0 0 0 +3 0 1 1 +4 0 1 0 +5 1 0 1 +6 1 0 0 +7 1 1 1 +8 1 1 0 +9 0 0.5 0.5 +10 1 0.5 0.5 +11 0.5 0 0.5 +12 0.5 1 0.5 +13 0.5 0.5 0 +14 0.5 0.5 1 +$EndNodes +$Elements +68 +1 15 2 0 1 1 +2 15 2 0 2 2 +3 15 2 0 3 3 +4 15 2 0 4 4 +5 15 2 0 5 5 +6 15 2 0 6 6 +7 15 2 0 7 7 +8 15 2 0 8 8 +9 1 2 0 1 2 1 +10 1 2 0 2 1 3 +11 1 2 0 3 4 3 +12 1 2 0 4 2 4 +13 1 2 0 5 6 5 +14 1 2 0 6 5 7 +15 1 2 0 7 8 7 +16 1 2 0 8 6 8 +17 1 2 0 9 2 6 +18 1 2 0 10 1 5 +19 1 2 0 11 4 8 +20 1 2 0 12 3 7 +21 2 2 0 1 2 1 9 +22 2 2 0 1 1 3 9 +23 2 2 0 1 4 2 9 +24 2 2 0 1 3 4 9 +25 2 2 0 2 6 10 5 +26 2 2 0 2 5 10 7 +27 2 2 0 2 8 10 6 +28 2 2 0 2 7 10 8 +29 2 2 0 3 1 2 11 +30 2 2 0 3 5 1 11 +31 2 2 0 3 2 6 11 +32 2 2 0 3 6 5 11 +33 2 2 0 4 3 12 4 +34 2 2 0 4 7 12 3 +35 2 2 0 4 4 12 8 +36 2 2 0 4 8 12 7 +37 2 2 0 5 2 4 13 +38 2 2 0 5 6 2 13 +39 2 2 0 5 4 8 13 +40 2 2 0 5 8 6 13 +41 2 2 0 6 1 14 3 +42 2 2 0 6 5 14 1 +43 2 2 0 6 3 14 7 +44 2 2 0 6 7 14 5 +45 4 2 0 1 13 10 9 11 +46 4 2 0 1 14 9 12 10 +47 4 2 0 1 10 9 12 13 +48 4 2 0 1 14 9 10 11 +49 4 2 0 1 12 7 10 8 +50 4 2 0 1 9 12 4 3 +51 4 2 0 1 5 1 14 11 +52 4 2 0 1 14 3 9 1 +53 4 2 0 1 12 3 14 7 +54 4 2 0 1 13 8 10 6 +55 4 2 0 1 7 5 14 10 +56 4 2 0 1 9 2 4 13 +57 4 2 0 1 1 2 9 11 +58 4 2 0 1 12 4 8 13 +59 4 2 0 1 6 13 2 11 +60 4 2 0 1 6 5 10 11 +61 4 2 0 1 9 14 1 11 +62 4 2 0 1 3 9 12 14 +63 4 2 0 1 14 12 7 10 +64 4 2 0 1 4 12 9 13 +65 4 2 0 1 8 10 12 13 +66 4 2 0 1 10 13 6 11 +67 4 2 0 1 10 5 14 11 +68 4 2 0 1 9 2 13 11 +$EndElements