slate kernels with dS_h wrong for 1 layer meshes

[wence-]

no kernel should be executed, but one is.

from firedrake import *
mesh = UnitTriangleMesh()
mesh = ExtrudedMesh(mesh, 1)
DG = FunctionSpace(mesh, "DG", 0) 
n = FacetNormal(mesh) 
u = TestFunction(DG) 

x, y, z = SpatialCoordinate(mesh) 
form = jump(z*u)*dS_h 

A = assemble(Tensor(form)).dat.data 
ref = assemble(form).dat.data   

A => array([0.5])
ref => array([0.])

[sv2518]

Do you know what the problem is?

[wence-]

There's an off-by-one in the checking of which kernels to execute for the interior horizontal kernels.

[sv2518]

Can I just do the following?

diff --git a/firedrake/slate/slac/compiler.py b/firedrake/slate/slac/compiler.py
index 064acaec..93a7fc72 100644
--- a/firedrake/slate/slac/compiler.py
+++ b/firedrake/slate/slac/compiler.py
@@ -497,7 +497,8 @@ def tensor_assembly_calls(builder):
         bottom = ast.Block(int_top + ext_btm, open_scope=True)
         top = ast.Block(int_btm + ext_top, open_scope=True)
         rest = ast.Block(int_btm + int_top, open_scope=True)
-        statements.append(ast.If(ast.Eq(builder.mesh_layer_sym, 0),
+        check = 1 if (int_top or int_btm) else 1
+        statements.append(ast.If(ast.Eq(builder.mesh_layer_sym, check),
                                  (bottom, ast.If(eq_layer, (top, rest)))))

[wence-]

I think the problem is that the logic is "am I on the bottom?", OK, then execute the "int_top" kernel. But it should be something like ("am I on the bottom? OK, am I on the top? -> do nothing; ...).

FWIW, both branches your condition result in check = 1.

[sv2518]

Ah oops that was only for a quick test for what happens when it's 1 and I use dS_v. I actually had 0 in the else.

[wence-]

I think it would be easier to split these conditionals into a sequence of ifs (rather the current if else setup).

if (current_layer == 0) 
    ext_bottom;
if (top_layer)
   ext_top;
if (current_layer == 0 && !top_layer)
   int_top;
if (top_layer && current_layer != 0)
   int_btm;
if (current_layer > 0 && !top_layer)
   int_top;
   int_btm

Or something?

[wence-]

The problem with the current set of conditions is that it merges kernels that come from exterior face terms and those that come from interior face terms. Whereas one needs to handle them separately, I think. For example, your proposed fix would fix my bug, but would break if you had a ds_b (rather than dS_h integral) I think.

[sv2518]

I cannot get a failing example for dS_v. Do I need to change the direction of extrusion for that?

[wence-]

dS_v will not fail, since it goes down a different codepath.

[sv2518]

Sorry, I don't think I understand your suggestion. What is top_layer? Current_layer we have available over a symbol via builder.mesh_layer_sym, but I do not think we have top_layer? Did you mean int_top or ext_top here https://github.com/firedrakeproject/firedrake/blob/2a8aea3ffe3d0868fa218cb5c70f555d239e23e9/firedrake/slate/slac/compiler.py#L492?

[wence-]

What is top_layer? Current_layer we have available over a symbol via builder.mesh_layer_sym, but I do not think we have top_layer?

Oh, it was pseudo-code for "conditional checking that we're on the top layer".

[sv2518]

Ah, how do I check that? current_layer == layer_count?

[wence-]

Or maybe current_layer == layer_count - 1 ?

[sv2518]

If we determine the top_layer over the current_layer, that makes your proposed conditions redundant though.

[wence-]

OK, so let's step back, what are the cases we actually need to distinguish, and what are the conditions we need to make sure the right kernels are executed? I think I vaguely sketched something that might be right, but clearly it is not fully correct. Perhaps one needs a sequence of if-elseif-...

[sv2518]

I counter proposed a fix here 3704a83. It's the way I do it in TSSLAC.