diff --git a/.pylintrc-local.yml b/.pylintrc-local.yml
index be1948583b86a81c56b102f4c159e1282e5d31dd..08f36d4c97f33574ae9c40bded4b8fba226de85b 100644
--- a/.pylintrc-local.yml
+++ b/.pylintrc-local.yml
@@ -9,3 +9,6 @@
- pml.py
- poisson.py
- second_order.py
+- arg: ignored-modules
+ val:
+ - sympy
diff --git a/doc/utils.rst b/doc/utils.rst
index 2235db8f518489069dc52b20683cbe60e838d457..a6f068b6f8ef033d4641f72b2810ce4889ac50e3 100644
--- a/doc/utils.rst
+++ b/doc/utils.rst
@@ -10,3 +10,8 @@ Array contexts
--------------
.. automodule:: grudge.array_context
+
+Miscellaneous tools
+-------------------
+
+.. automodule:: grudge.tools
diff --git a/grudge/models/wave.py b/grudge/models/wave.py
index 1cb1cdb4a2813f35c621adbea8b814f5d97879ce..209575ba3258f5cfcf983dc49b13f880ef2c8748 100644
--- a/grudge/models/wave.py
+++ b/grudge/models/wave.py
@@ -182,6 +182,9 @@ class WeakWaveOperator(HyperbolicOperator):
def max_characteristic_velocity(self, actx, t=None, fields=None):
return abs(self.c)
+ def estimate_rk4_timestep(self, actx, dcoll, **kwargs):
+ # FIXME: Sketchy, empirically determined fudge factor
+ return 0.38 * super().estimate_rk4_timestep(actx, dcoll, **kwargs)
# }}}
diff --git a/grudge/tools.py b/grudge/tools.py
index 25018b530ccf1748f20c4a04191bcacf8e622799..9723f9249414c7b268c064d938944100991db2aa 100644
--- a/grudge/tools.py
+++ b/grudge/tools.py
@@ -1,4 +1,6 @@
-"""Miscellaneous helper facilities."""
+"""
+.. autofunction:: build_jacobian
+"""
__copyright__ = "Copyright (C) 2007 Andreas Kloeckner"
@@ -199,3 +201,38 @@ class OrderedSet(MutableSet):
return set(self) == set(other)
# }}}
+
+
+def build_jacobian(actx, f, base_state, stepsize):
+ """Returns a Jacobian matrix of *f* determined by a one-sided finite
+ difference approximation with *stepsize*.
+
+ :param f: a callable with a single argument, any array or
+ :class:`arraycontext.ArrayContainer` supported by *actx*.
+ :param base_state: The point at which the Jacobian is to be
+ calculated. May be any array or :class:`arraycontext.ArrayContainer`
+ supported by *actx*.
+ :returns: a two-dimensional :class:`numpy.ndarray`
+ """
+ from arraycontext import flatten, unflatten
+ flat_base_state = flatten(base_state, actx)
+
+ n, = flat_base_state.shape
+
+ mat = np.empty((n, n), dtype=flat_base_state.dtype)
+
+ f_base = f(base_state)
+
+ for i in range(n):
+ unit_i_flat = np.zeros(n, dtype=mat.dtype)
+ unit_i_flat[i] = stepsize
+
+ f_unit_i = f(f_base + unflatten(
+ base_state, actx.from_numpy(unit_i_flat), actx))
+
+ mat[:, i] = actx.to_numpy(flatten((f_unit_i - f_base) / stepsize, actx))
+
+ return mat
+
+
+# vim: foldmethod=marker
diff --git a/requirements.txt b/requirements.txt
index abf6d39fa3d227e72056061b4fa08efa6d9044a5..2107e5aeb28297c113c8e4b92fcf7c11439df48a 100644
--- a/requirements.txt
+++ b/requirements.txt
@@ -15,3 +15,6 @@ git+https://github.com/inducer/pyvisfile.git#egg=pyvisfile
git+https://github.com/inducer/pymetis.git#egg=pymetis
git+https://github.com/illinois-ceesd/logpyle.git#egg=logpyle
git+https://github.com/inducer/pytato.git#egg=pytato
+
+# for test_wave_dt_estimate
+sympy
diff --git a/test/test_dt_utils.py b/test/test_dt_utils.py
index 573e7bad6c6c390000b929846ccfaf6fcdd23a9b..af9cbb6b4b88384d568ad6aa3158914cb7577fdf 100644
--- a/test/test_dt_utils.py
+++ b/test/test_dt_utils.py
@@ -45,6 +45,7 @@ import logging
logger = logging.getLogger(__name__)
+from meshmode import _acf # noqa: F401
@pytest.mark.parametrize("name", ["interval", "box2d", "box3d"])
@@ -123,6 +124,74 @@ def test_non_geometric_factors(actx_factory, name):
assert all(factors <= upper_bounds)
+@pytest.mark.parametrize("dim", [1, 2])
+@pytest.mark.parametrize("degree", [2, 4])
+def test_wave_dt_estimate(actx_factory, dim, degree, visualize=False):
+ actx = actx_factory()
+
+ 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)
+ assert mesh.dim == dim
+
+ dcoll = DiscretizationCollection(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))
+
+ from pytools.obj_array import make_obj_array
+ fields = make_obj_array([dcoll.zeros(actx) for i in range(dim+1)])
+
+ from grudge.tools import build_jacobian
+ mat = build_jacobian(actx, rhs, fields, 1)
+
+ import numpy.linalg as la
+ eigvals = la.eigvals(mat)
+
+ assert (eigvals.real <= 1e-12).all()
+
+ from leap.rk import stability_function, RK4MethodBuilder
+ import sympy as sp
+ stab_func = sp.lambdify(*stability_function(
+ RK4MethodBuilder.a_explicit,
+ RK4MethodBuilder.output_coeffs))
+
+ dt_est = actx.to_numpy(wave_op.estimate_rk4_timestep(actx, dcoll))
+
+ if visualize:
+ re, im = np.mgrid[-4:1:30j, -5:5:30j]
+ sf_grid = np.abs(stab_func(re+1j*im))
+
+ 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.grid()
+ plt.show()
+
+ thresh = 1+1e-8
+ max_stab = np.max(np.abs(stab_func(dt_est*eigvals)))
+ assert max_stab < thresh, max_stab
+
+ dt_factors = 2**np.linspace(0, 4, 40)[1:]
+ stable_dt_factors = [
+ dt_factor
+ for dt_factor in dt_factors
+ if np.max(np.abs(stab_func(dt_factor*dt_est*eigvals))) < thresh]
+
+ if stable_dt_factors:
+ 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
+
+
# You can test individual routines by typing
# $ python test_grudge.py 'test_routine()'
diff --git a/test/test_mpi_communication.py b/test/test_mpi_communication.py
index 82a6cc234afdf1e0e0cccd8c182ca1274362c9ae..24e3d63805be26d573e0292c06109c9a2ddc6f95 100644
--- a/test/test_mpi_communication.py
+++ b/test/test_mpi_communication.py
@@ -243,9 +243,8 @@ def _test_mpi_wave_op_entrypoint(actx, visualize=False):
[dcoll.zeros(actx) for i in range(dcoll.dim)]
)
- # FIXME: Sketchy, empirically determined fudge factor
dt = actx.to_numpy(
- 0.45 * wave_op.estimate_rk4_timestep(actx, dcoll, fields=fields))
+ wave_op.estimate_rk4_timestep(actx, dcoll, fields=fields))
wave_op.check_bc_coverage(local_mesh)