From d2c80ca8dbce9ab1b71fdefa6b0737c90da59837 Mon Sep 17 00:00:00 2001
From: Andreas Kloeckner <inform@tiker.net>
Date: Sun, 19 Feb 2017 13:05:40 -0600
Subject: [PATCH] Make weak-form operators work
---
grudge/execution.py | 19 +++++++--------
grudge/models/advection.py | 46 +++++++++++++++++++-----------------
grudge/symbolic/operators.py | 15 ++++++++----
test/test_grudge.py | 18 ++++++++++----
4 files changed, 57 insertions(+), 41 deletions(-)
diff --git a/grudge/execution.py b/grudge/execution.py
index 2fa47b91..ca6b7dd9 100644
--- a/grudge/execution.py
+++ b/grudge/execution.py
@@ -343,10 +343,6 @@ class ExecutionMapper(mappers.Evaluator,
def exec_diff_batch_assign(self, insn):
field = self.rec(insn.field)
repr_op = insn.operators[0]
- if not isinstance(repr_op, sym.RefDiffOperator):
- # FIXME
- raise NotImplementedError()
-
# FIXME: There's no real reason why differentiation is special,
# execution-wise.
# This should be unified with map_elementwise_linear, which should
@@ -376,21 +372,24 @@ class ExecutionMapper(mappers.Evaluator,
noperators = len(insn.operators)
result = discr.empty(
- dtype=field.dtype, extra_dims=(noperators,)).with_queue(self.queue)
+ queue=self.queue,
+ dtype=field.dtype, extra_dims=(noperators,),
+ allocator=self.bound_op.allocator)
for grp in discr.groups:
if grp.nelements == 0:
continue
- # FIXME: Should make matrices on device and cache them
- diff_matrices = grp.diff_matrices()
- diff_matrices_ary = np.empty((
+ matrices = repr_op.matrices(grp)
+
+ # FIXME: Should transfer matrices to device and cache them
+ matrices_ary = np.empty((
noperators, grp.nunit_nodes, grp.nunit_nodes))
for i, op in enumerate(insn.operators):
- diff_matrices_ary[i] = diff_matrices[op.rst_axis]
+ matrices_ary[i] = matrices[op.rst_axis]
knl()(self.queue,
- diff_mat=diff_matrices_ary,
+ diff_mat=matrices_ary,
result=grp.view(result), vec=grp.view(field))
return [(name, result[i]) for i, name in enumerate(insn.names)], []
diff --git a/grudge/models/advection.py b/grudge/models/advection.py
index da920cc0..8fda51a0 100644
--- a/grudge/models/advection.py
+++ b/grudge/models/advection.py
@@ -88,41 +88,43 @@ class StrongAdvectionOperator(AdvectionOperatorBase):
self.flux(pair))
return (
- -self.v.dot(sym.nabla(self.ambient_dim)*u)
+ - self.v.dot(sym.nabla(self.ambient_dim)*u)
+ sym.InverseMassOperator()(
sym.FaceMassOperator()(
flux(sym.int_tpair(u))
- + flux(sym.bv_tpair(sym.BTAG_ALL, u, self.inflow_u)))))
+ + flux(sym.bv_tpair(sym.BTAG_ALL, u, self.inflow_u))
+
+ # FIXME: Add back support for inflow/outflow tags
+ #+ flux(sym.bv_tpair(self.inflow_tag, u, bc_in))
+ #+ flux(sym.bv_tpair(self.outflow_tag, u, bc_out))
+ )))
class WeakAdvectionOperator(AdvectionOperatorBase):
- def flux(self):
- return self.weak_flux()
+ def flux(self, u):
+ return self.weak_flux(u)
def sym_operator(self):
- from grudge.symbolic import (
- Field,
- BoundaryPair,
- get_flux_operator,
- make_stiffness_t,
- InverseMassOperator,
- RestrictToBoundary)
-
- u = Field("u")
+ u = sym.var("u")
# boundary conditions -------------------------------------------------
- bc_in = Field("bc_in")
- bc_out = RestrictToBoundary(self.outflow_tag)*u
+ bc_in = self.inflow_u
+ # bc_out = sym.interp("vol", self.outflow_tag)(u)
- stiff_t = make_stiffness_t(self.dimensions)
- m_inv = InverseMassOperator()
+ def flux(pair):
+ return sym.interp(pair.dd, "all_faces")(
+ self.flux(pair))
- flux_op = get_flux_operator(self.flux())
+ return sym.InverseMassOperator()(
+ np.dot(
+ self.v, sym.stiffness_t(self.ambient_dim)*u)
+ - sym.FaceMassOperator()(
+ flux(sym.int_tpair(u))
+ + flux(sym.bv_tpair(sym.BTAG_ALL, u, bc_in))
- return m_inv(np.dot(self.v, stiff_t*u) - (
- flux_op(u)
- + flux_op(BoundaryPair(u, bc_in, self.inflow_tag))
- + flux_op(BoundaryPair(u, bc_out, self.outflow_tag))
+ # FIXME: Add back support for inflow/outflow tags
+ #+ flux(sym.bv_tpair(self.inflow_tag, u, bc_in))
+ #+ flux(sym.bv_tpair(self.outflow_tag, u, bc_out))
))
# }}}
diff --git a/grudge/symbolic/operators.py b/grudge/symbolic/operators.py
index 69759298..dc2e4fa1 100644
--- a/grudge/symbolic/operators.py
+++ b/grudge/symbolic/operators.py
@@ -214,10 +214,20 @@ class RefDiffOperatorBase(ElementwiseLinearOperator):
class RefDiffOperator(RefDiffOperatorBase):
mapper_method = intern("map_ref_diff")
+ @staticmethod
+ def matrices(element_group):
+ return element_group.diff_matrices()
+
class RefStiffnessTOperator(RefDiffOperatorBase):
mapper_method = intern("map_ref_stiffness_t")
+ @staticmethod
+ def matrices(element_group):
+ assert element_group.is_orthogonal_basis()
+ mmat = element_group.mass_matrix()
+ return [dmat.T.dot(mmat.T) for dmat in element_group.diff_matrices()]
+
# }}}
# }}}
@@ -349,10 +359,7 @@ class RefMassOperatorBase(ElementwiseLinearOperator):
class RefMassOperator(RefMassOperatorBase):
@staticmethod
def matrix(element_group):
- import modepy as mp
- return mp.mass_matrix(
- element_group.basis(),
- element_group.unit_nodes)
+ return element_group.mass_matrix()
mapper_method = intern("map_ref_mass")
diff --git a/test/test_grudge.py b/test/test_grudge.py
index b4b6db45..33725674 100644
--- a/test/test_grudge.py
+++ b/test/test_grudge.py
@@ -208,10 +208,11 @@ def test_2d_gauss_theorem(ctx_factory):
("rect2", [4, 8]),
("rect3", [4, 6]),
])
+@pytest.mark.parametrize("op_type", ["strong", "weak"])
@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):
+def test_convergence_advec(ctx_factory, mesh_name, mesh_pars, op_type, flux_type,
+ order, visualize=False):
"""Test whether 2D advection actually converges"""
cl_ctx = cl.create_some_context()
@@ -268,13 +269,20 @@ def test_convergence_advec(ctx_factory, mesh_name, mesh_pars, flux_type, order,
-v.dot(x)/norm_v
+ sym.var("t", sym.DD_SCALAR)*norm_v)
- from grudge.models.advection import StrongAdvectionOperator
+ from grudge.models.advection import (
+ StrongAdvectionOperator, WeakAdvectionOperator)
discr = Discretization(cl_ctx, mesh, order=order)
- op = StrongAdvectionOperator(v,
+ op_class = {
+ "strong": StrongAdvectionOperator,
+ "weak": WeakAdvectionOperator,
+ }[op_type]
+ op = op_class(v,
inflow_u=u_analytic(sym.nodes(dim, sym.BTAG_ALL)),
flux_type=flux_type)
bound_op = bind(discr, op.sym_operator())
+ #print(bound_op)
+ #1/0
u = bind(discr, u_analytic(sym.nodes(dim)))(queue, t=0)
@@ -309,7 +317,7 @@ def test_convergence_advec(ctx_factory, mesh_name, mesh_pars, flux_type, order,
last_u = event.state_component
if visualize:
- vis.write_vtk_file("fld-%04d.vtu" % step,
+ vis.write_vtk_file("fld-%s-%04d.vtu" % (mesh_par, step),
[("u", event.state_component)])
error_l2 = bind(discr,
--
GitLab