diff --git a/examples/advection/surface.py b/examples/advection/surface.py index 56bdfa3d2e98382e7138146393464fc6804ebdfa..9dcdc93ec550acdb0e50d210ce470c897bb04d24 100644 --- a/examples/advection/surface.py +++ b/examples/advection/surface.py @@ -203,7 +203,7 @@ def main(ctx_factory, dim=2, order=4, use_quad=False, visualize=False): # {{{ time stepping - dt = adv_operator.estimate_rk4_timestep(dcoll, fields=u0, dt_scaling=8/9) + dt = adv_operator.estimate_rk4_timestep(dcoll, fields=u0) nsteps = int(final_time // dt) + 1 logger.info("dt: %.5e", dt) diff --git a/examples/wave/var-propagation-speed.py b/examples/wave/var-propagation-speed.py index a0237d7bd6b814d7fb0391ddefb99387404b703f..0b8d51b9bdfa0388612366c1d95cedcc931eedb2 100644 --- a/examples/wave/var-propagation-speed.py +++ b/examples/wave/var-propagation-speed.py @@ -102,7 +102,8 @@ def main(ctx_factory, dim=2, order=4, visualize=False): def rhs(t, w): return wave_op.operator(t, w) - dt = wave_op.estimate_rk4_timestep(dcoll, fields=fields) + dt_scaling_const = 2/3 + dt = dt_scaling_const * wave_op.estimate_rk4_timestep(dcoll, fields=fields) dt_stepper = set_up_rk4("w", dt, fields, rhs) diff --git a/examples/wave/wave-min-mpi.py b/examples/wave/wave-min-mpi.py index 25f0174dc57f31d058cd224f8f4a5da80e9c3bd2..a3e75c7b525628bcaa6bd861a1009ec9cd79d985 100644 --- a/examples/wave/wave-min-mpi.py +++ b/examples/wave/wave-min-mpi.py @@ -114,7 +114,8 @@ def main(ctx_factory, dim=2, order=4, visualize=False): [dcoll.zeros(actx) for i in range(dcoll.dim)] ) - dt = wave_op.estimate_rk4_timestep(dcoll, fields=fields) + dt_scaling_const = 2/3 + dt = dt_scaling_const * wave_op.estimate_rk4_timestep(dcoll, fields=fields) wave_op.check_bc_coverage(local_mesh) diff --git a/examples/wave/wave-min.py b/examples/wave/wave-min.py index 7f19a2287e4f27e7e75cd646bbb63f22e093c262..7b18d3dfaec68ab38af936c45f395bd10ad29a75 100644 --- a/examples/wave/wave-min.py +++ b/examples/wave/wave-min.py @@ -95,7 +95,8 @@ def main(ctx_factory, dim=2, order=4, visualize=False): [dcoll.zeros(actx) for i in range(dcoll.dim)] ) - dt = wave_op.estimate_rk4_timestep(dcoll, fields=fields) + dt_scaling_const = 2/3 + dt = dt_scaling_const * wave_op.estimate_rk4_timestep(dcoll, fields=fields) wave_op.check_bc_coverage(mesh) diff --git a/examples/wave/wave-op-mpi.py b/examples/wave/wave-op-mpi.py index 4c0d01bc31ea57b0f65f72a3c5717ff06197ad5b..3fb49d7d1b65eb8d9caa5ac6a5d68c79625c5a94 100644 --- a/examples/wave/wave-op-mpi.py +++ b/examples/wave/wave-op-mpi.py @@ -114,14 +114,14 @@ def rk4_step(y, t, h, f): return y + h/6*(k1 + 2*k2 + 2*k3 + k4) -def estimate_rk4_timestep(dcoll, c, scaling=None): +def estimate_rk4_timestep(dcoll, c): from grudge.dt_utils import (dt_non_geometric_factor, dt_geometric_factors) # FIXME: We should make the nodal reductions respect MPI. # See: https://github.com/inducer/grudge/issues/112 and # https://github.com/inducer/grudge/issues/113 - dt_factor = (dt_non_geometric_factor(dcoll, scaling=scaling) + dt_factor = (dt_non_geometric_factor(dcoll) * op.nodal_min(dcoll, "vol", dt_geometric_factors(dcoll))) mpi_comm = dcoll.mpi_communicator @@ -192,7 +192,8 @@ def main(ctx_factory, dim=2, order=3, visualize=False): ) c = 1 - dt = estimate_rk4_timestep(dcoll, c, scaling=0.45) + dt_scaling_const = 0.45 + dt = dt_scaling_const * estimate_rk4_timestep(dcoll, c) vis = make_visualizer(dcoll) diff --git a/examples/wave/wave-op-var-velocity.py b/examples/wave/wave-op-var-velocity.py index 1f5c35b3b60d86ff40214a0ae4432e15a0058a85..e48c8d7141285b40d0a0935c032449e7bbe5a2a4 100644 --- a/examples/wave/wave-op-var-velocity.py +++ b/examples/wave/wave-op-var-velocity.py @@ -129,11 +129,11 @@ def rk4_step(y, t, h, f): return y + h/6*(k1 + 2*k2 + 2*k3 + k4) -def estimate_rk4_timestep(dcoll, c, scaling=None): +def estimate_rk4_timestep(dcoll, c): from grudge.dt_utils import (dt_non_geometric_factor, dt_geometric_factors) - dt_factor = (dt_non_geometric_factor(dcoll, scaling=scaling) + dt_factor = (dt_non_geometric_factor(dcoll) * op.nodal_min(dcoll, "vol", dt_geometric_factors(dcoll))) return dt_factor * (1 / c) @@ -189,7 +189,8 @@ def main(ctx_factory, dim=2, order=3, visualize=False): # bounded above by 1 c = 0.2 + 0.8*bump(actx, dcoll, center=np.zeros(3), width=0.5) - dt = estimate_rk4_timestep(dcoll, c=1, scaling=0.5) + dt_scaling_const = 0.5 + dt = dt_scaling_const * estimate_rk4_timestep(dcoll, c=1) fields = flat_obj_array( bump(actx, dcoll, ), diff --git a/examples/wave/wave-op.py b/examples/wave/wave-op.py index 9c39f289904688131e131822ffb1c8507b338156..c0494d9f0196da5cb7d951a6950f82250a2e9aa6 100644 --- a/examples/wave/wave-op.py +++ b/examples/wave/wave-op.py @@ -112,11 +112,11 @@ def rk4_step(y, t, h, f): return y + h/6*(k1 + 2*k2 + 2*k3 + k4) -def estimate_rk4_timestep(dcoll, c, scaling=None): +def estimate_rk4_timestep(dcoll, c): from grudge.dt_utils import (dt_non_geometric_factor, dt_geometric_factors) - dt_factor = (dt_non_geometric_factor(dcoll, scaling=scaling) + dt_factor = (dt_non_geometric_factor(dcoll) * op.nodal_min(dcoll, "vol", dt_geometric_factors(dcoll))) return dt_factor * (1 / c) @@ -165,7 +165,8 @@ def main(ctx_factory, dim=2, order=3, visualize=False): ) c = 1 - dt = estimate_rk4_timestep(dcoll, c, scaling=0.45) + dt_scaling_const = 0.45 + dt = dt_scaling_const * estimate_rk4_timestep(dcoll, c) vis = make_visualizer(dcoll) diff --git a/grudge/dt_utils.py b/grudge/dt_utils.py index 4fdba9e7db5ce6d9362bba27d298f615e98fec86..14c955595fc5ca41ddcd707066e0238865b273e5 100644 --- a/grudge/dt_utils.py +++ b/grudge/dt_utils.py @@ -1,7 +1,7 @@ """Helper functions for estimating stable time steps for RKDG methods. .. autofunction:: dt_non_geometric_factor -.. autofunction:: dt_geometric_factor +.. autofunction:: dt_geometric_factors """ __copyright__ = """ @@ -44,8 +44,7 @@ from pytools import memoize_on_first_arg @memoize_on_first_arg -def dt_non_geometric_factor( - dcoll: DiscretizationCollection, scaling=None, dd=None) -> float: +def dt_non_geometric_factor(dcoll: DiscretizationCollection, dd=None) -> float: r"""Computes the non-geometric scale factor following [Hesthaven_2008]_, section 6.4: @@ -56,8 +55,6 @@ def dt_non_geometric_factor( where :math:`\Delta r_i` denotes the distance between two distinct nodes on the reference element. - :arg scaling: a :class:`float` denoting the scaling factor. By default, - the constant is set to 2/3. :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 :class:`float` denoting the minimum node distance on the @@ -66,9 +63,6 @@ def dt_non_geometric_factor( if dd is None: dd = DD_VOLUME - if scaling is None: - scaling = 2/3 - discr = dcoll.discr_from_dd(dd) min_delta_rs = [] for grp in discr.groups: @@ -93,7 +87,7 @@ def dt_non_geometric_factor( ) # Return minimum over all element groups in the discretization - return scaling * min(min_delta_rs) + return min(min_delta_rs) @memoize_on_first_arg @@ -132,12 +126,14 @@ def dt_geometric_factors( "Geometric factors are only implemented for simplex element groups" ) - cell_vols = op.elementwise_integral( - dcoll, dd, volm_discr.zeros(actx) + 1.0 + cell_vols = abs( + op.elementwise_integral( + dcoll, dd, volm_discr.zeros(actx) + 1.0 + ) ) if dcoll.dim == 1: - return op.nodal_min(dcoll, dd, cell_vols) + return cell_vols dd_face = DOFDesc("all_faces", dd.discretization_tag) face_discr = dcoll.discr_from_dd(dd_face) @@ -146,8 +142,10 @@ def dt_geometric_factors( # take the sum over the averaged face areas on each face. # NOTE: The face areas are the *same* at each face nodal location. # This assumes there are the *same* number of face nodes on each face. - surface_areas = op.elementwise_integral( - dcoll, dd_face, face_discr.zeros(actx) + 1.0 + surface_areas = abs( + op.elementwise_integral( + dcoll, dd_face, face_discr.zeros(actx) + 1.0 + ) ) surface_areas = DOFArray( actx, @@ -162,7 +160,7 @@ def dt_geometric_factors( for vgrp, afgrp, face_ae_i in zip(volm_discr.groups, face_discr.groups, - actx.np.fabs(surface_areas)) + surface_areas) ) ) @@ -177,4 +175,3 @@ def dt_geometric_factors( for cv_i, sae_i in zip(cell_vols, surface_areas) ) ) -