From 8f553da5e62aa9bc90da82f6c6d26a5b6d51f64e Mon Sep 17 00:00:00 2001
From: Andreas Kloeckner <inform@tiker.net>
Date: Mon, 24 May 2021 17:10:33 -0500
Subject: [PATCH] Drop model tests from old symbolic tests
---
test/test_grudge_sym_old.py | 300 ------------------------------------
1 file changed, 300 deletions(-)
diff --git a/test/test_grudge_sym_old.py b/test/test_grudge_sym_old.py
index f246a895..ed47028f 100644
--- a/test/test_grudge_sym_old.py
+++ b/test/test_grudge_sym_old.py
@@ -21,7 +21,6 @@ THE SOFTWARE.
"""
import numpy as np
-import numpy.linalg as la
from meshmode import _acf # noqa: F401
from meshmode.dof_array import flatten, thaw
@@ -651,305 +650,6 @@ def test_surface_divergence_theorem(actx_factory, mesh_name, visualize=False):
# }}}
-# {{{ models: advection
-
-@pytest.mark.parametrize(("mesh_name", "mesh_pars"), [
- ("segment", [8, 16, 32]),
- ("disk", [0.1, 0.05]),
- ("rect2", [4, 8]),
- ("rect3", [4, 6]),
- ("warped2", [4, 8]),
- ])
-@pytest.mark.parametrize("op_type", ["strong", "weak"])
-@pytest.mark.parametrize("flux_type", ["central"])
-@pytest.mark.parametrize("order", [3, 4, 5])
-# test: 'test_convergence_advec(cl._csc, "disk", [0.1, 0.05], "strong", "upwind", 3)'
-def test_convergence_advec(actx_factory, mesh_name, mesh_pars, op_type, flux_type,
- order, visualize=False):
- """Test whether 2D advection actually converges"""
-
- actx = actx_factory()
-
- from pytools.convergence import EOCRecorder
- eoc_rec = EOCRecorder()
-
- for mesh_par in mesh_pars:
- if mesh_name == "segment":
- mesh = mgen.generate_box_mesh(
- [np.linspace(-1.0, 1.0, mesh_par)],
- order=order)
-
- dim = 1
- dt_factor = 1.0
- elif mesh_name == "disk":
- pytest.importorskip("meshpy")
-
- from meshpy.geometry import make_circle, GeometryBuilder
- from meshpy.triangle import MeshInfo, build
-
- geob = GeometryBuilder()
- geob.add_geometry(*make_circle(1))
- mesh_info = MeshInfo()
- geob.set(mesh_info)
-
- mesh_info = build(mesh_info, max_volume=mesh_par)
-
- from meshmode.mesh.io import from_meshpy
- mesh = from_meshpy(mesh_info, order=1)
- dim = 2
- dt_factor = 4
- elif mesh_name.startswith("rect"):
- dim = int(mesh_name[-1:])
- mesh = mgen.generate_regular_rect_mesh(a=(-0.5,)*dim, b=(0.5,)*dim,
- nelements_per_axis=(mesh_par,)*dim, order=4)
-
- if dim == 2:
- dt_factor = 4
- elif dim == 3:
- dt_factor = 2
- else:
- raise ValueError("dt_factor not known for %dd" % dim)
- elif mesh_name.startswith("warped"):
- dim = int(mesh_name[-1:])
- mesh = mgen.generate_warped_rect_mesh(dim, order=order,
- nelements_side=mesh_par)
-
- if dim == 2:
- dt_factor = 4
- elif dim == 3:
- dt_factor = 2
- else:
- raise ValueError("dt_factor not known for %dd" % dim)
- else:
- raise ValueError("invalid mesh name: " + mesh_name)
-
- v = np.array([0.27, 0.31, 0.1])[:dim]
- norm_v = la.norm(v)
-
- def f(x):
- return sym.sin(10*x)
-
- def u_analytic(x):
- return f(
- -v.dot(x)/norm_v
- + sym.var("t", dof_desc.DD_SCALAR)*norm_v)
-
- from grudge.models.advection import (
- StrongAdvectionOperator, WeakAdvectionOperator)
- from meshmode.mesh import BTAG_ALL
-
- discr = DiscretizationCollection(actx, mesh, order=order)
- op_class = {
- "strong": StrongAdvectionOperator,
- "weak": WeakAdvectionOperator,
- }[op_type]
- op = op_class(v,
- inflow_u=u_analytic(sym.nodes(dim, BTAG_ALL)),
- flux_type=flux_type)
-
- bound_op = bind(discr, op.sym_operator())
-
- u = bind(discr, u_analytic(sym.nodes(dim)))(actx, t=0)
-
- def rhs(t, u):
- return bound_op(t=t, u=u)
-
- if dim == 3:
- final_time = 0.1
- else:
- final_time = 0.2
-
- h_max = bind(discr, sym.h_max_from_volume(discr.ambient_dim))(actx)
- dt = dt_factor * h_max/order**2
- nsteps = (final_time // dt) + 1
- dt = final_time/nsteps + 1e-15
-
- from grudge.shortcuts import set_up_rk4
- dt_stepper = set_up_rk4("u", dt, u, rhs)
-
- last_u = None
-
- from grudge.shortcuts import make_visualizer
- vis = make_visualizer(discr, vis_order=order)
-
- step = 0
-
- for event in dt_stepper.run(t_end=final_time):
- if isinstance(event, dt_stepper.StateComputed):
- step += 1
- logger.debug("[%04d] t = %.5f", step, event.t)
-
- last_t = event.t
- last_u = event.state_component
-
- if visualize:
- vis.write_vtk_file("fld-%s-%04d.vtu" % (mesh_par, step),
- [("u", event.state_component)])
-
- error_l2 = bind(discr,
- sym.norm(2, sym.var("u")-u_analytic(sym.nodes(dim))))(
- t=last_t, u=last_u)
- logger.info("h_max %.5e error %.5e", h_max, error_l2)
- eoc_rec.add_data_point(h_max, error_l2)
-
- logger.info("\n%s", eoc_rec.pretty_print(
- abscissa_label="h",
- error_label="L2 Error"))
-
- if mesh_name.startswith("warped"):
- # NOTE: curvilinear meshes are hard
- assert eoc_rec.order_estimate() > order - 0.5
- else:
- assert eoc_rec.order_estimate() > order
-
-# }}}
-
-
-# {{{ models: maxwell
-
-@pytest.mark.parametrize("order", [3, 4, 5])
-def test_convergence_maxwell(actx_factory, order):
- """Test whether 3D Maxwell's actually converges"""
-
- actx = actx_factory()
-
- from pytools.convergence import EOCRecorder
- eoc_rec = EOCRecorder()
-
- dims = 3
- ns = [4, 6, 8]
- for n in ns:
- mesh = mgen.generate_regular_rect_mesh(
- a=(0.0,)*dims,
- b=(1.0,)*dims,
- nelements_per_axis=(n,)*dims)
-
- discr = DiscretizationCollection(actx, mesh, order=order)
-
- epsilon = 1
- mu = 1
-
- from grudge.models.em import get_rectangular_cavity_mode
- sym_mode = get_rectangular_cavity_mode(1, (1, 2, 2))
-
- analytic_sol = bind(discr, sym_mode)
- fields = analytic_sol(actx, t=0, epsilon=epsilon, mu=mu)
-
- from grudge.models.em import MaxwellOperator
- op = MaxwellOperator(epsilon, mu, flux_type=0.5, dimensions=dims)
- op.check_bc_coverage(mesh)
- bound_op = bind(discr, op.sym_operator())
-
- def rhs(t, w):
- return bound_op(t=t, w=w)
-
- dt = 0.002
- final_t = dt * 5
- nsteps = int(final_t/dt)
-
- from grudge.shortcuts import set_up_rk4
- dt_stepper = set_up_rk4("w", dt, fields, rhs)
-
- logger.info("dt %.5e nsteps %5d", dt, nsteps)
-
- norm = bind(discr, sym.norm(2, sym.var("u")))
-
- step = 0
- for event in dt_stepper.run(t_end=final_t):
- if isinstance(event, dt_stepper.StateComputed):
- assert event.component_id == "w"
- esc = event.state_component
-
- step += 1
- logger.debug("[%04d] t = %.5e", step, event.t)
-
- sol = analytic_sol(actx, mu=mu, epsilon=epsilon, t=step * dt)
- vals = [norm(u=(esc[i] - sol[i])) / norm(u=sol[i]) for i in range(5)] # noqa E501
- total_error = sum(vals)
- eoc_rec.add_data_point(1.0/n, total_error)
-
- logger.info("\n%s", eoc_rec.pretty_print(
- abscissa_label="h",
- error_label="L2 Error"))
-
- assert eoc_rec.order_estimate() > order
-
-# }}}
-
-
-# {{{ models: variable coefficient advection oversampling
-
-@pytest.mark.parametrize("order", [2, 3, 4])
-def test_improvement_quadrature(actx_factory, order):
- """Test whether quadrature improves things and converges"""
- from grudge.models.advection import VariableCoefficientAdvectionOperator
- from pytools.convergence import EOCRecorder
- from meshmode.discretization.poly_element import QuadratureSimplexGroupFactory
-
- actx = actx_factory()
-
- dims = 2
- sym_nds = sym.nodes(dims)
- advec_v = flat_obj_array(-1*sym_nds[1], sym_nds[0])
-
- flux = "upwind"
- op = VariableCoefficientAdvectionOperator(advec_v, 0, flux_type=flux)
-
- def gaussian_mode():
- source_width = 0.1
- sym_x = sym.nodes(2)
- return sym.exp(-np.dot(sym_x, sym_x) / source_width**2)
-
- def conv_test(descr, use_quad):
- logger.info("-" * 75)
- logger.info(descr)
- logger.info("-" * 75)
- eoc_rec = EOCRecorder()
-
- ns = [20, 25]
- for n in ns:
- mesh = mgen.generate_regular_rect_mesh(
- a=(-0.5,)*dims,
- b=(0.5,)*dims,
- nelements_per_axis=(n,)*dims,
- order=order)
-
- if use_quad:
- discr_tag_to_group_factory = {
- "product": QuadratureSimplexGroupFactory(order=4*order)
- }
- else:
- discr_tag_to_group_factory = {"product": None}
-
- discr = DiscretizationCollection(
- actx, mesh, order=order,
- discr_tag_to_group_factory=discr_tag_to_group_factory
- )
-
- bound_op = bind(discr, op.sym_operator())
- fields = bind(discr, gaussian_mode())(actx, t=0)
- norm = bind(discr, sym.norm(2, sym.var("u")))
-
- esc = bound_op(u=fields)
- total_error = norm(u=esc)
- eoc_rec.add_data_point(1.0/n, total_error)
-
- logger.info("\n%s", eoc_rec.pretty_print(
- abscissa_label="h",
- error_label="L2 Error"))
-
- return eoc_rec.order_estimate(), np.array([x[1] for x in eoc_rec.history])
-
- eoc, errs = conv_test("no quadrature", False)
- q_eoc, q_errs = conv_test("with quadrature", True)
-
- assert q_eoc > eoc
- assert (q_errs < errs).all()
- assert q_eoc > order - 0.1
-
-# }}}
-
-
# {{{ operator collector determinism
def test_op_collector_order_determinism():
--
GitLab