diff --git a/grudge/models/wave.py b/grudge/models/wave.py index b22ff159..9f07797d 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 aa6ebc09..04375a37 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