From 168574c1c9f7dfdc46c409ebbc48efa4d66f84fb Mon Sep 17 00:00:00 2001
From: Alexandru Fikl <alexfikl@gmail.com>
Date: Sun, 3 May 2020 21:26:39 -0500
Subject: [PATCH] use h_max in advection convergence test
---
grudge/function_registry.py | 1 +
grudge/symbolic/operators.py | 7 +++++--
grudge/symbolic/primitives.py | 1 -
test/test_grudge.py | 15 +++++++--------
4 files changed, 13 insertions(+), 11 deletions(-)
diff --git a/grudge/function_registry.py b/grudge/function_registry.py
index 2f56eba4..4d521560 100644
--- a/grudge/function_registry.py
+++ b/grudge/function_registry.py
@@ -136,6 +136,7 @@ class CElementwiseBinaryFunction(Function):
return func(arg0, arg1)
from pymbolic.primitives import Variable
+
@memoize_in(self, "map_call_knl_%s" % func_name)
def knl():
i = Variable("i")
diff --git a/grudge/symbolic/operators.py b/grudge/symbolic/operators.py
index 2352b189..1a573058 100644
--- a/grudge/symbolic/operators.py
+++ b/grudge/symbolic/operators.py
@@ -695,17 +695,20 @@ def norm(p, arg, dd=None):
raise ValueError("unsupported value of p")
-def h_max(dd=None):
+def h_max(ambient_dim, dim=None, dd=None):
import grudge.symbolic.primitives as prim
if dd is None:
dd = prim.DD_VOLUME
dd = prim.as_dofdesc(dd)
+ if dim is None:
+ dim = ambient_dim
+
return NodalMax(dd_in=dd)(
ElementwiseSumOperator(dd)(
MassOperator(dd_in=dd)(prim.Ones(dd))
)
- )
+ )**(1.0/dim)
# }}}
diff --git a/grudge/symbolic/primitives.py b/grudge/symbolic/primitives.py
index 8ff71cc4..aa00e8e6 100644
--- a/grudge/symbolic/primitives.py
+++ b/grudge/symbolic/primitives.py
@@ -654,7 +654,6 @@ def mv_normal(dd, ambient_dim, dim=None):
def normal(dd, ambient_dim, dim=None, quadrature_tag=None):
return mv_normal(dd, ambient_dim, dim).as_vector()
-
# }}}
diff --git a/test/test_grudge.py b/test/test_grudge.py
index 65cd9d89..e27599d9 100644
--- a/test/test_grudge.py
+++ b/test/test_grudge.py
@@ -232,7 +232,6 @@ def test_convergence_advec(ctx_factory, mesh_name, mesh_pars, op_type, flux_type
[np.linspace(-1.0, 1.0, mesh_par)],
order=order)
- h = 2.0 / mesh_par
dim = 1
dt_factor = 1.0
elif mesh_name == "disk":
@@ -250,7 +249,6 @@ def test_convergence_advec(ctx_factory, mesh_name, mesh_pars, op_type, flux_type
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"):
@@ -259,7 +257,6 @@ def test_convergence_advec(ctx_factory, mesh_name, mesh_pars, op_type, flux_type
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:
@@ -306,7 +303,8 @@ def test_convergence_advec(ctx_factory, mesh_name, mesh_pars, op_type, flux_type
else:
final_time = 0.2
- dt = dt_factor * h/order**2
+ h_max = bind(discr, sym.h_max(discr.ambient_dim))(queue)
+ dt = dt_factor * h_max/order**2
nsteps = (final_time // dt) + 1
dt = final_time/nsteps + 1e-15
@@ -335,11 +333,12 @@ def test_convergence_advec(ctx_factory, mesh_name, mesh_pars, op_type, flux_type
error_l2 = bind(discr,
sym.norm(2, sym.var("u")-u_analytic(sym.nodes(dim))))(
queue, t=last_t, u=last_u)
- print(h, error_l2)
- eoc_rec.add_data_point(h, error_l2)
+ print(h_max, error_l2)
+ eoc_rec.add_data_point(h_max, error_l2)
- print(eoc_rec.pretty_print(abscissa_label="h",
- error_label="L2 Error"))
+ print(eoc_rec.pretty_print(
+ abscissa_label="h",
+ error_label="L2 Error"))
assert eoc_rec.order_estimate() > order
--
GitLab