Skip to content

Commit

Permalink
Include warpy mesh in dt test.
Browse files Browse the repository at this point in the history
  • Loading branch information
MTCam committed Jul 28, 2024
1 parent 24e17ae commit 7b8f37b
Show file tree
Hide file tree
Showing 2 changed files with 32 additions and 19 deletions.
2 changes: 1 addition & 1 deletion grudge/models/wave.py
Original file line number Diff line number Diff line change
Expand Up @@ -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)

Expand Down
49 changes: 31 additions & 18 deletions test/test_dt_utils.py
Original file line number Diff line number Diff line change
Expand Up @@ -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"])
Expand Down Expand Up @@ -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)])
Expand Down Expand Up @@ -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()

Expand All @@ -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
Expand Down

0 comments on commit 7b8f37b

Please sign in to comment.