Make weak-form operators work

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)], []

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))
                     ))
 
 # }}}

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")
 

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,