Incorrect exterior vertical boundary integral with variable layers extrusion

[stephankramer]

The following code

from firedrake import *
from math import ceil

L = 100
H1 = 2.
H2 = 42.
dy = 5.0
ny = round(L/dy)
dz = 2.0

# create mesh
mesh1d = IntervalMesh(ny, L)
layers = []
cell = 0
yr = 0
for i in range(ny):
    yr += dy  # y of rhs of column (assumed to be the higher one)
    height = H1 + yr/L * (H2-H1)
    ncells = ceil(height/dz)
    layers.append([0, ncells])
    cell += ncells

mesh = ExtrudedMesh(mesh1d, layers, layer_height=dz)
y = mesh.coordinates.dat.data_ro[:,0]
z = mesh.coordinates.dat.data_ro[:,1]
mesh.coordinates.dat.data[:,1] = np.minimum(H1 + y/L * (H2-H1), z)

print("Length of bottom", assemble(Constant(1.0)*ds_b(domain=mesh)))
print("Length of top", assemble(Constant(1.0)*ds_t(domain=mesh)))
print("Length of lhs", assemble(Constant(1.0)*ds_v(1, domain=mesh)))
print("Length of rhs", assemble(Constant(1.0)*ds_v(2, domain=mesh)))

sets up an extruded domain for the ocean under an iceshelf (2d vertical slice, flat bottom, left boundary of 2m, right boundary of 42m, with sloping top boundary). The first three boundary integrals (bottom, top and lhs) are calculated correctly. The last one comes out incorrectly (giving 4 instead of the expected 42) where it seems that the rhs boundary only has 2 vertical facets (the same as in the first column on the left) instead of 22 vertical facets. This is confirmed by looking at something like solve(v*u*dx-v*Constant(1.0)*ds_v(2, domain=mesh)==0, u) with a DG Function u and TestFunction v, where the solution only appears in the bottom right corner. If I look at boundary nodes using DirichletBC it does come up with the correct nodes along the rhs boundary.

[wence-]

Try this patch in pyop2?

diff --git a/pyop2/codegen/builder.py b/pyop2/codegen/builder.py
index 81811ea8..cbb2b111 100644
--- a/pyop2/codegen/builder.py
+++ b/pyop2/codegen/builder.py
@@ -702,8 +702,6 @@ class WrapperBuilder(object):
     def _layer_index(self):
         if self.constant_layers:
             return FixedIndex(0)
-        if self.subset:
-            return self._loop_index
         else:
             return self.loop_index

[stephankramer]

That works - magic!

[wence-]

Do you mind preparing a test?

[stephankramer]

Thanks! Sure, will do