From 811b3c09d16e6445d38f27bfc39bf5f8cf605bd0 Mon Sep 17 00:00:00 2001 From: Mike Campbell Date: Thu, 13 Jun 2024 15:29:08 -0500 Subject: [PATCH 1/4] Add TPE option, reorder term calc in euler to expose issue with face normals in lazy --- examples/euler/acoustic_pulse.py | 28 +++++++++++++++++----------- grudge/models/euler.py | 12 ++++++------ 2 files changed, 23 insertions(+), 17 deletions(-) diff --git a/examples/euler/acoustic_pulse.py b/examples/euler/acoustic_pulse.py index 779062910..d5bbc0375 100644 --- a/examples/euler/acoustic_pulse.py +++ b/examples/euler/acoustic_pulse.py @@ -37,6 +37,7 @@ EulerOperator, InviscidWallBC ) +from meshmode.mesh import TensorProductElementGroup from grudge.shortcuts import rk4_step from meshmode.mesh import BTAG_ALL @@ -112,7 +113,8 @@ def run_acoustic_pulse(actx, final_time=1, resolution=16, overintegration=False, - visualize=False): + visualize=False, + tpe=False): # eos-related parameters gamma = 1.4 @@ -124,16 +126,18 @@ def run_acoustic_pulse(actx, dim = 2 box_ll = -0.5 box_ur = 0.5 + group_cls = TensorProductElementGroup if tpe else None mesh = generate_regular_rect_mesh( a=(box_ll,)*dim, b=(box_ur,)*dim, - nelements_per_axis=(resolution,)*dim) + nelements_per_axis=(resolution,)*dim, + group_cls=group_cls) from grudge import DiscretizationCollection from grudge.dof_desc import DISCR_TAG_BASE, DISCR_TAG_QUAD - from meshmode.discretization.poly_element import \ - (default_simplex_group_factory, - QuadratureSimplexGroupFactory) + from meshmode.discretization.poly_element import ( + InterpolatoryEdgeClusteredGroupFactory, + QuadratureGroupFactory) exp_name = f"fld-acoustic-pulse-N{order}-K{resolution}" if overintegration: @@ -145,9 +149,8 @@ def run_acoustic_pulse(actx, dcoll = DiscretizationCollection( actx, mesh, discr_tag_to_group_factory={ - DISCR_TAG_BASE: default_simplex_group_factory( - base_dim=mesh.dim, order=order), - DISCR_TAG_QUAD: QuadratureSimplexGroupFactory(2*order) + DISCR_TAG_BASE: InterpolatoryEdgeClusteredGroupFactory(order), + DISCR_TAG_QUAD: QuadratureGroupFactory(2*order) } ) @@ -212,7 +215,8 @@ def rhs(t, q): def main(ctx_factory, order=3, final_time=1, resolution=16, - overintegration=False, visualize=False, lazy=False): + overintegration=False, visualize=False, lazy=False, + tpe=False): cl_ctx = ctx_factory() queue = cl.CommandQueue(cl_ctx) @@ -234,7 +238,7 @@ def main(ctx_factory, order=3, final_time=1, resolution=16, resolution=resolution, overintegration=overintegration, final_time=final_time, - visualize=visualize + visualize=visualize, tpe=tpe ) @@ -251,6 +255,8 @@ def main(ctx_factory, order=3, final_time=1, resolution=16, help="write out vtk output") parser.add_argument("--lazy", action="store_true", help="switch to a lazy computation mode") + parser.add_argument("--tpe", action="store_true", + help="use tensor product elements") args = parser.parse_args() logging.basicConfig(level=logging.INFO) @@ -260,4 +266,4 @@ def main(ctx_factory, order=3, final_time=1, resolution=16, resolution=args.resolution, overintegration=args.oi, visualize=args.visualize, - lazy=args.lazy) + lazy=args.lazy, tpe=args.tpe) diff --git a/grudge/models/euler.py b/grudge/models/euler.py index 1b6eb569c..25d26e780 100644 --- a/grudge/models/euler.py +++ b/grudge/models/euler.py @@ -322,12 +322,6 @@ def operator(self, t, q): def interp_to_quad(u): return op.project(dcoll, "vol", dq, u) - # Compute volume fluxes - volume_fluxes = op.weak_local_div( - dcoll, dq, - interp_to_quad(euler_volume_flux(dcoll, q, gamma=gamma)) - ) - # Compute interior interface fluxes interface_fluxes = ( sum( @@ -357,6 +351,12 @@ def interp_to_quad(u): ) interface_fluxes = interface_fluxes + bc_fluxes + # Compute volume fluxes + volume_fluxes = op.weak_local_div( + dcoll, dq, + interp_to_quad(euler_volume_flux(dcoll, q, gamma=gamma)) + ) + return op.inverse_mass( dcoll, volume_fluxes - op.face_mass(dcoll, df, interface_fluxes) From 5b563de6fc84d7fed2d543944bdd89410aeb9d64 Mon Sep 17 00:00:00 2001 From: Mike Campbell Date: Mon, 17 Jun 2024 17:39:59 -0500 Subject: [PATCH 2/4] Use geoderiv_connection only for Simplices --- grudge/geometry/metrics.py | 7 ++++++- 1 file changed, 6 insertions(+), 1 deletion(-) diff --git a/grudge/geometry/metrics.py b/grudge/geometry/metrics.py index 260f219b3..3e8177cca 100644 --- a/grudge/geometry/metrics.py +++ b/grudge/geometry/metrics.py @@ -86,6 +86,11 @@ register_multivector_as_array_container() +def _has_geoderiv_connection(grp): + from modepy.shapes import Simplex + return grp.is_affine and issubclass(grp._modepy_shape_cls, Simplex) + + def _geometry_to_quad_if_requested( dcoll, inner_dd, dd, vec, _use_geoderiv_connection): @@ -105,7 +110,7 @@ def to_quad(vec): return DOFArray( vec.array_context, tuple( - geoderiv_vec_i if megrp.is_affine else all_quad_vec_i + geoderiv_vec_i if _has_geoderiv_connection(megrp) else all_quad_vec_i for megrp, geoderiv_vec_i, all_quad_vec_i in zip( dcoll.discr_from_dd(inner_dd).mesh.groups, dcoll._base_to_geoderiv_connection(inner_dd)(vec), From 6df78a7c57a2e4bf0e3aaf9fb89d19762649a94c Mon Sep 17 00:00:00 2001 From: Mike Campbell Date: Mon, 17 Jun 2024 17:40:54 -0500 Subject: [PATCH 3/4] Allow generation of TPE mesh type --- test/mesh_data.py | 16 +++++++++++++++- 1 file changed, 15 insertions(+), 1 deletion(-) diff --git a/test/mesh_data.py b/test/mesh_data.py index c71950bb9..6658f495e 100644 --- a/test/mesh_data.py +++ b/test/mesh_data.py @@ -4,6 +4,7 @@ from meshmode.mesh.io import read_gmsh import numpy as np import meshmode.mesh.generation as mgen +from meshmode.mesh import TensorProductElementGroup class MeshBuilder(ABC): @@ -111,6 +112,7 @@ def get_mesh(self, resolution, mesh_order=4): class _BoxMeshBuilderBase(MeshBuilder): resolutions = [4, 8, 16] mesh_order = 1 + group_cls = None a = (-0.5, -0.5, -0.5) b = (+0.5, +0.5, +0.5) @@ -122,20 +124,32 @@ def get_mesh(self, resolution, mesh_order=4): return mgen.generate_regular_rect_mesh( a=self.a, b=self.b, nelements_per_axis=resolution, - order=mesh_order) + order=mesh_order, group_cls=self.group_cls) class BoxMeshBuilder1D(_BoxMeshBuilderBase): ambient_dim = 1 + def __init__(self, tpe=False): + if tpe: + self.group_cls = TensorProductElementGroup + class BoxMeshBuilder2D(_BoxMeshBuilderBase): ambient_dim = 2 + def __init__(self, tpe=False): + if tpe: + self.group_cls = TensorProductElementGroup + class BoxMeshBuilder3D(_BoxMeshBuilderBase): ambient_dim = 2 + def __init__(self, tpe=False): + if tpe: + self.group_cls = TensorProductElementGroup + class WarpedRectMeshBuilder(MeshBuilder): resolutions = [4, 6, 8] From 00b6ca45efc2b24daa6caf4190addc75945ed34b Mon Sep 17 00:00:00 2001 From: Mike Campbell Date: Mon, 17 Jun 2024 17:41:26 -0500 Subject: [PATCH 4/4] Extend some tests to hit TPEs --- test/test_grudge.py | 48 +++++++++++++++++++++++++++++++-------------- 1 file changed, 33 insertions(+), 15 deletions(-) diff --git a/test/test_grudge.py b/test/test_grudge.py index ff42eb63f..c1432f37d 100644 --- a/test/test_grudge.py +++ b/test/test_grudge.py @@ -61,10 +61,10 @@ # {{{ mass operator trig integration +@pytest.mark.parametrize("tpe", [True, False]) @pytest.mark.parametrize("ambient_dim", [1, 2, 3]) -@pytest.mark.parametrize("discr_tag", [dof_desc.DISCR_TAG_BASE, - dof_desc.DISCR_TAG_QUAD]) -def test_mass_mat_trig(actx_factory, ambient_dim, discr_tag): +@pytest.mark.parametrize("use_overint", [False, True]) +def test_mass_mat_trig(actx_factory, tpe, ambient_dim, use_overint): """Check the integral of some trig functions on an interval using the mass matrix. """ @@ -75,22 +75,26 @@ def test_mass_mat_trig(actx_factory, ambient_dim, discr_tag): a = -4.0 * np.pi b = +9.0 * np.pi + true_integral = 13*np.pi/2 * (b - a)**(ambient_dim - 1) + discr_tag = dof_desc.DISCR_TAG_QUAD if use_overint else dof_desc.DISCR_TAG_BASE - from meshmode.discretization.poly_element import QuadratureSimplexGroupFactory dd_quad = dof_desc.DOFDesc(dof_desc.DTAG_VOLUME_ALL, discr_tag) + discr_order = order if discr_tag is dof_desc.DISCR_TAG_BASE: discr_tag_to_group_factory = {} else: + discr_order = None discr_tag_to_group_factory = { - discr_tag: QuadratureSimplexGroupFactory(order=2*order) + dof_desc.DISCR_TAG_BASE: InterpolatoryEdgeClusteredGroupFactory(order), + dof_desc.DISCR_TAG_QUAD: QuadratureGroupFactory(order=2*order) } mesh = mgen.generate_regular_rect_mesh( a=(a,)*ambient_dim, b=(b,)*ambient_dim, nelements_per_axis=(nel_1d,)*ambient_dim, order=1) dcoll = DiscretizationCollection( - actx, mesh, order=order, + actx, mesh, order=discr_order, discr_tag_to_group_factory=discr_tag_to_group_factory ) @@ -160,10 +164,11 @@ def _spheroid_surface_area(radius, aspect_ratio): return 2.0 * np.pi * radius**2 * (1 + (c/a)**2 / e * np.arctanh(e)) +@pytest.mark.parametrize("tpe", [True, False]) @pytest.mark.parametrize("name", [ "2-1-ellipse", "spheroid", "box2d", "box3d" ]) -def test_mass_surface_area(actx_factory, name): +def test_mass_surface_area(actx_factory, tpe, name): actx = actx_factory() # {{{ cases @@ -171,16 +176,20 @@ def test_mass_surface_area(actx_factory, name): order = 4 if name == "2-1-ellipse": + if tpe: + pytest.skip() builder = mesh_data.EllipseMeshBuilder(radius=3.1, aspect_ratio=2.0) surface_area = _ellipse_surface_area(builder.radius, builder.aspect_ratio) elif name == "spheroid": + if tpe: + pytest.skip() builder = mesh_data.SpheroidMeshBuilder() surface_area = _spheroid_surface_area(builder.radius, builder.aspect_ratio) elif name == "box2d": - builder = mesh_data.BoxMeshBuilder2D() + builder = mesh_data.BoxMeshBuilder2D(tpe) surface_area = 1.0 elif name == "box3d": - builder = mesh_data.BoxMeshBuilder3D() + builder = mesh_data.BoxMeshBuilder3D(tpe) surface_area = 1.0 else: raise ValueError("unknown geometry name: %s" % name) @@ -976,11 +985,11 @@ def rhs(t, w): # {{{ models: variable coefficient advection oversampling @pytest.mark.parametrize("order", [2, 3, 4]) -def test_improvement_quadrature(actx_factory, order): +@pytest.mark.parametrize("tpe", [False, True]) +def test_improvement_quadrature(actx_factory, order, tpe): """Test whether quadrature improves things and converges""" from grudge.models.advection import VariableCoefficientAdvectionOperator from pytools.convergence import EOCRecorder - from meshmode.discretization.poly_element import QuadratureSimplexGroupFactory from meshmode.mesh import BTAG_ALL actx = actx_factory() @@ -1002,23 +1011,30 @@ def conv_test(descr, use_quad): else: qtag = None + group_cls = TensorProductElementGroup if tpe else None + qfac = 2 if tpe else 4 ns = [20, 25] + discr_order = order for n in ns: mesh = mgen.generate_regular_rect_mesh( a=(-0.5,)*dims, b=(0.5,)*dims, nelements_per_axis=(n,)*dims, - order=order) + order=order, group_cls=group_cls) if use_quad: discr_tag_to_group_factory = { - qtag: QuadratureSimplexGroupFactory(order=4*order) + dof_desc.DISCR_TAG_BASE: + InterpolatoryEdgeClusteredGroupFactory(order), + dof_desc.DISCR_TAG_QUAD: + QuadratureGroupFactory(order=qfac*order) } + discr_order = None else: discr_tag_to_group_factory = {} dcoll = DiscretizationCollection( - actx, mesh, order=order, + actx, mesh, order=discr_order, discr_tag_to_group_factory=discr_tag_to_group_factory ) @@ -1050,9 +1066,11 @@ def zero_inflow(dtag, t=0): eoc, errs = conv_test("no quadrature", False) q_eoc, q_errs = conv_test("with quadrature", True) - assert q_eoc > eoc assert (q_errs < errs).all() assert q_eoc > order - 0.1 + # Fails for all tensor-product element types + assert q_eoc > eoc + # }}}