From 6fc72f523747e001a2ebcb5a224b8fa4b68658e9 Mon Sep 17 00:00:00 2001
From: Andreas Kloeckner <inform@tiker.net>
Date: Sun, 10 Jan 2016 18:47:34 -0600
Subject: [PATCH] Add advection convergence test
---
test/test_grudge.py | 125 +++++++++++++++++++++++++++++++++++++++++++-
1 file changed, 124 insertions(+), 1 deletion(-)
diff --git a/test/test_grudge.py b/test/test_grudge.py
index ee1828a5..18d4aba4 100644
--- a/test/test_grudge.py
+++ b/test/test_grudge.py
@@ -164,6 +164,8 @@ def test_tri_diff_mat(ctx_factory, dim, order=4):
def test_2d_gauss_theorem(ctx_factory):
"""Verify Gauss's theorem explicitly on a mesh"""
+ pytest.importorskip("meshpy")
+
from meshpy.geometry import make_circle, GeometryBuilder
from meshpy.triangle import MeshInfo, build
@@ -201,8 +203,129 @@ def test_2d_gauss_theorem(ctx_factory):
assert abs(gauss_err.get()) < 5e-15
+@pytest.mark.parametrize(("mesh_name", "mesh_pars"), [
+ ("disk", [0.1, 0.05]),
+ ("rect2", [4, 8]),
+ ("rect3", [4, 6]),
+ ])
+@pytest.mark.parametrize("flux_type", ["upwind"])
+@pytest.mark.parametrize("order", [3, 4, 5])
+def test_convergence_advec(ctx_factory, mesh_name, mesh_pars, flux_type, order,
+ visualize=False):
+ """Test whether 2D advection actually converges"""
+
+ cl_ctx = cl.create_some_context()
+ queue = cl.CommandQueue(cl_ctx)
+
+ from pytools.convergence import EOCRecorder
+ eoc_rec = EOCRecorder()
+
+ for mesh_par in mesh_pars:
+ if 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)
+ h = np.sqrt(mesh_par)
+ dim = 2
+ dt_factor = 4
+
+ elif mesh_name.startswith("rect"):
+ dim = int(mesh_name[4:])
+ from meshmode.mesh.generation import generate_regular_rect_mesh
+ mesh = generate_regular_rect_mesh(a=(-0.5,)*dim, b=(0.5,)*dim,
+ n=(mesh_par,)*dim, order=4)
+
+ h = 1/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", sym.DD_SCALAR)*norm_v)
+
+ from grudge.models.advection import StrongAdvectionOperator
+ discr = Discretization(cl_ctx, mesh, order=order)
+ op = StrongAdvectionOperator(v,
+ inflow_u=u_analytic(sym.nodes(dim, sym.BTAG_ALL)),
+ flux_type=flux_type)
+
+ bound_op = bind(discr, op.sym_operator())
+
+ u = bind(discr, u_analytic(sym.nodes(dim)))(queue, t=0)
+
+ def rhs(t, u):
+ return bound_op(queue, t=t, u=u)
+
+ if dim == 3:
+ final_time = 0.1
+ else:
+ final_time = 0.2
+
+ dt = dt_factor * h/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
+ print(event.t)
+
+ last_t = event.t
+ last_u = event.state_component
+
+ if visualize:
+ vis.write_vtk_file("fld-%04d.vtu" % step,
+ [("u", event.state_component)])
+
+ error_l2 = bind(discr,
+ sym.norm(2, sym.var("u")-u_analytic(sym.nodes(dim))))(
+ queue, t=last_t, u=last_u).get()
+ print(h, error_l2)
+ eoc_rec.add_data_point(h, error_l2)
+
+ print(eoc_rec.pretty_print(abscissa_label="h",
+ error_label="L2 Error"))
+
+ assert eoc_rec.order_estimate() > order
+
+
# You can test individual routines by typing
-# $ python test_layer_pot.py 'test_routine()'
+# $ python test_grudge.py 'test_routine()'
if __name__ == "__main__":
import sys
--
GitLab