diff --git a/examples/gas-in-box.py b/examples/gas-in-box.py index b1fa426b2..bfac79294 100644 --- a/examples/gas-in-box.py +++ b/examples/gas-in-box.py @@ -47,7 +47,7 @@ from mirgecom.utils import force_evaluation from mirgecom.io import make_init_message -from mirgecom.integrators import rk4_step +from mirgecom.integrators import rk4_step, euler_step from mirgecom.steppers import advance_state from mirgecom.boundary import ( AdiabaticSlipBoundary, @@ -102,8 +102,8 @@ def main(actx_class, use_esdg=False, use_tpe=False, mech_name="uiuc_7sp", transport_type=0, use_av=0, use_limiter=False, order=1, nscale=1, npassive_species=0, map_mesh=False, - rotation_angle=0, add_pulse=False, - mesh_filename=None): + rotation_angle=0, add_pulse=False, nsteps=20, + mesh_filename=None, euler_timestepping=False): """Drive the example.""" if casename is None: casename = "gas-in-box" @@ -132,11 +132,11 @@ def main(actx_class, use_esdg=False, use_tpe=False, from leap.rk import RK4MethodBuilder timestepper = RK4MethodBuilder("state") else: - timestepper = rk4_step - n_steps = 20 + timestepper = euler_step if euler_timestepping else rk4_step + current_cfl = 1.0 current_dt = 1e-6 - t_final = current_dt * n_steps + t_final = current_dt * nsteps current_t = 0 constant_cfl = False temperature_tolerance = 1e-2 @@ -193,6 +193,7 @@ def main(actx_class, use_esdg=False, use_tpe=False, local_mesh, global_nelements = distribute_mesh(comm, generate_mesh) local_nelements = local_mesh.nelements + if dim is None: dim = local_mesh.ambient_dim @@ -221,8 +222,7 @@ def add_wonk(x: np.ndarray) -> np.ndarray: theta = rotation_angle/180.0 * np.pi local_mesh = rotate_mesh_around_axis(local_mesh, theta=theta) - dcoll = create_discretization_collection(actx, local_mesh, order=order, - tensor_product_elements=use_tpe) + dcoll = create_discretization_collection(actx, local_mesh, order=order) nodes = actx.thaw(dcoll.nodes()) ones = dcoll.zeros(actx) + 1. @@ -751,6 +751,10 @@ def my_rhs(t, stepper_state): help="use leap timestepper") parser.add_argument("--esdg", action="store_true", help="use entropy-stable dg for inviscid terms.") + parser.add_argument("--euler-timestepping", action="store_true", + help="use euler timestepping") + parser.add_argument("--nsteps", type=int, default=20, + help="number of timesteps to take") parser.add_argument("--numpy", action="store_true", help="use numpy-based eager actx.") parser.add_argument("-d", "--dimension", type=int, choices=[1, 2, 3], @@ -818,7 +822,7 @@ def my_rhs(t, stepper_state): main(actx_class, use_esdg=args.esdg, dim=args.dimension, use_overintegration=args.overintegration or args.esdg, - use_leap=args.leap, use_tpe=args.tpe, + use_leap=args.leap, use_tpe=args.tpe, nsteps=args.nsteps, casename=args.casename, rst_filename=rst_filename, periodic_mesh=args.periodic, use_mixture=args.mixture, multiple_boundaries=args.boundaries, @@ -828,6 +832,6 @@ def my_rhs(t, stepper_state): use_navierstokes=args.navierstokes, npassive_species=args.species, nscale=args.weak_scale, mech_name=args.mechanism_name, map_mesh=args.wonky, rotation_angle=args.rotate, add_pulse=args.pulse, - mesh_filename=args.meshfile) + mesh_filename=args.meshfile, euler_timestepping=args.euler_timestepping) # vim: foldmethod=marker diff --git a/mirgecom/discretization.py b/mirgecom/discretization.py index 11b48b318..2d4458678 100644 --- a/mirgecom/discretization.py +++ b/mirgecom/discretization.py @@ -49,35 +49,26 @@ def create_discretization_collection(actx, volume_meshes, order, *, DeprecationWarning, stacklevel=2) if tensor_product_elements: - warn("Overintegration is not supported for tensor product elements.") + warn( + "tensor_product_elements argument is not needed and will vanish soon.", + DeprecationWarning, stacklevel=2) from grudge.dof_desc import DISCR_TAG_BASE, DISCR_TAG_QUAD, DISCR_TAG_MODAL from grudge.discretization import make_discretization_collection from meshmode.discretization.poly_element import ( - QuadratureSimplexGroupFactory, - PolynomialRecursiveNodesGroupFactory, - LegendreGaussLobattoTensorProductGroupFactory as Lgl, + InterpolatoryEdgeClusteredGroupFactory, + QuadratureGroupFactory, ModalGroupFactory ) if quadrature_order < 0: - quadrature_order = 2*order+1 + quadrature_order = 2*order + 1 - if tensor_product_elements: - return make_discretization_collection( - actx, volume_meshes, - discr_tag_to_group_factory={ - DISCR_TAG_BASE: Lgl(order), - DISCR_TAG_MODAL: ModalGroupFactory(order) - } - ) - else: - return make_discretization_collection( - actx, volume_meshes, - discr_tag_to_group_factory={ - DISCR_TAG_BASE: PolynomialRecursiveNodesGroupFactory(order=order, - family="lgl"), - DISCR_TAG_QUAD: QuadratureSimplexGroupFactory(quadrature_order), - DISCR_TAG_MODAL: ModalGroupFactory(order) - } - ) + return make_discretization_collection( + actx, volume_meshes, + discr_tag_to_group_factory={ + DISCR_TAG_BASE: InterpolatoryEdgeClusteredGroupFactory(order), + DISCR_TAG_MODAL: ModalGroupFactory(order), + DISCR_TAG_QUAD: QuadratureGroupFactory(quadrature_order) + } + ) diff --git a/test/test_operators.py b/test/test_operators.py index f08a8fb41..04ee9e86d 100644 --- a/test/test_operators.py +++ b/test/test_operators.py @@ -75,17 +75,17 @@ def _coord_test_func(dim, order=1): def _trig_test_func(dim): """Make trig test function. - 1d: cos(2pi x) - 2d: sin(2pi x)cos(2pi y) - 3d: sin(2pi x)sin(2pi y)cos(2pi z) + 1d: cos(pi x/2) + 2d: sin(pi x/2)cos(pi y/2) + 3d: sin(pi x/2)sin(pi y/2)cos(pi z/2) """ sym_coords = prim.make_sym_vector("x", dim) sym_cos = pmbl.var("cos") sym_sin = pmbl.var("sin") - sym_result = sym_cos(2*np.pi*sym_coords[dim-1]) + sym_result = sym_cos(np.pi*sym_coords[dim-1]/2) for i in range(dim-1): - sym_result = sym_result * sym_sin(2*np.pi*sym_coords[i]) + sym_result = sym_result * sym_sin(np.pi*sym_coords[i]/2) return sym_result @@ -133,8 +133,22 @@ def central_flux_boundary(actx, dcoll, soln_func, dd_bdry): return op.project(dcoll, bnd_tpair.dd, dd_allfaces, flux_weak) -@pytest.mark.parametrize("tpe", [False, True]) -@pytest.mark.parametrize("dim", [1, 2, 3]) +@pytest.mark.parametrize(("dim", "mesh_name", "rot_axis", "wonk"), + [ + (1, "tet_box1", None, False), + (2, "tet_box2", None, False), + (3, "tet_box3", None, False), + (2, "hex_box2", None, False), + (3, "hex_box3", None, False), + (2, "tet_box2_rot", np.array([0, 0, 1]), False), + (3, "tet_box3_rot1", np.array([0, 0, 1]), False), + (3, "tet_box3_rot2", np.array([0, 1, 1]), False), + (3, "tet_box3_rot3", np.array([1, 1, 1]), False), + (2, "hex_box2_rot", np.array([0, 0, 1]), False), + (3, "hex_box3_rot1", np.array([0, 0, 1]), False), + (3, "hex_box3_rot2", np.array([0, 1, 1]), False), + (3, "hex_box3_rot3", np.array([1, 1, 1]), False), + ]) @pytest.mark.parametrize("order", [1, 2, 3]) @pytest.mark.parametrize("sym_test_func_factory", [ partial(_coord_test_func, order=0), @@ -144,7 +158,8 @@ def central_flux_boundary(actx, dcoll, soln_func, dd_bdry): _trig_test_func, _cv_test_func ]) -def test_grad_operator(actx_factory, tpe, dim, order, sym_test_func_factory): +def test_grad_operator(actx_factory, dim, mesh_name, rot_axis, wonk, + order, sym_test_func_factory): """Test the gradient operator for sanity. Check whether we get the right answers for gradients of analytic functions with @@ -157,36 +172,67 @@ def test_grad_operator(actx_factory, tpe, dim, order, sym_test_func_factory): """ import pyopencl as cl from grudge.array_context import PyOpenCLArrayContext + from meshmode.mesh.processing import rotate_mesh_around_axis + from grudge.dt_utils import h_max_from_volume + from mirgecom.simutil import componentwise_norms + from arraycontext import flatten + from mirgecom.operators import grad_operator + from meshmode.mesh.processing import map_mesh + + def add_wonk(x: np.ndarray) -> np.ndarray: + wonk_field = np.empty_like(x) + if len(x) >= 2: + wonk_field[0] = ( + 1.5*x[0] + np.cos(x[0]) + + 0.1*np.sin(10*x[1])) + wonk_field[1] = ( + 0.05*np.cos(10*x[0]) + + 1.3*x[1] + np.sin(x[1])) + else: + wonk_field[0] = 1.5*x[0] + np.cos(x[0]) + + if len(x) >= 3: + wonk_field[2] = x[2] + np.sin(x[0] / 2) / 2 + + return wonk_field + + tpe = mesh_name.startswith("hex_") + rotation_angle = 32.0 + theta = rotation_angle/180.0 * np.pi # This comes from array_context actx = actx_factory() if tpe: # TPE requires *grudge* array context, not array_context - if dim == 1: # TPE does not support dim=1 - pytest.skip() ctx = cl.create_some_context() queue = cl.CommandQueue(ctx) actx = PyOpenCLArrayContext(queue) sym_test_func = sym_test_func_factory(dim) - tol = 1e-10 if dim < 3 else 1e-9 + tol = 1e-10 if dim < 3 else 2e-9 + from pytools.convergence import EOCRecorder eoc = EOCRecorder() - for nfac in [1, 2, 4]: + for nfac in [2, 4, 8]: + + # Make non-uniform spacings + b = tuple((2 ** n) for n in range(dim)) - mesh = get_box_mesh(dim, a=0, b=1, n=nfac*3, tensor_product_elements=tpe) + mesh = get_box_mesh(dim, a=0, b=b, n=nfac*3, tensor_product_elements=tpe) + if wonk: + mesh = map_mesh(mesh, add_wonk) + if rot_axis is not None: + mesh = rotate_mesh_around_axis(mesh, theta=theta, axis=rot_axis) logger.info( f"Number of {dim}d elements: {mesh.nelements}" ) - dcoll = create_discretization_collection(actx, mesh, order=order, - tensor_product_elements=tpe) + dcoll = create_discretization_collection(actx, mesh, order=order) # compute max element size - from grudge.dt_utils import h_max_from_volume h_max = h_max_from_volume(dcoll) def sym_eval(expr, x_vec): @@ -208,14 +254,13 @@ def sym_eval(expr, x_vec): test_data = test_func(nodes) exact_grad = grad_test_func(nodes) - from mirgecom.simutil import componentwise_norms - - from arraycontext import flatten err_scale = max(flatten(componentwise_norms(dcoll, exact_grad, np.inf), actx)) - if err_scale <= 1e-16: + print(f"{err_scale=}") + if err_scale <= 1e-10: err_scale = 1 + print(f"Rescaling: {err_scale=}") print(f"{test_data=}") print(f"{exact_grad=}") @@ -225,7 +270,6 @@ def sym_eval(expr, x_vec): test_data_flux_bnd = _elbnd_flux(dcoll, int_flux, bnd_flux, test_data_int_tpair, boundaries) - from mirgecom.operators import grad_operator dd_vol = as_dofdesc("vol") dd_allfaces = as_dofdesc("all_faces") test_grad = grad_operator(dcoll, dd_vol, dd_allfaces, @@ -236,6 +280,7 @@ def sym_eval(expr, x_vec): max(flatten( componentwise_norms(dcoll, test_grad - exact_grad, np.inf), actx) / err_scale) + print(f"{actx.to_numpy(h_max)=}\n{actx.to_numpy(grad_err)=}") eoc.add_data_point(actx.to_numpy(h_max), actx.to_numpy(grad_err))