Skip to content
Snippets Groups Projects
Commit 8f553da5 authored by Andreas Klöckner's avatar Andreas Klöckner
Browse files

Drop model tests from old symbolic tests

parent ff5c3ead
No related branches found
No related tags found
No related merge requests found
Hide whitespace changes
Inline Side-by-side
test/test_grudge_sym_old.py
+ 0
− 300
View file @ 8f553da5
......@@ -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():
......
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment