Slope limiter boundary treatment

[tkarna]

Extending slope limiters to work correctly at boundaries

[tkarna]

Added first implementation where the exterior facets are taken into account when computing the allowed min/max values. A couple of remarks:

When determining the allowed min/max values, this implementation uses the original field on the boundary facets. This means that the solution on the boundary facets is not limited at all. This is probably not what we want. We should therefore get a p0 solution on the facets somehow. For triangles one could compute the facet average in pyop2 if we knew the correct facet node indices. For quads the facet node average is not exactly the p0 solution but it might be good enough.

Second, If I try to use measure ds_t or ds_b in firedrake par_loop I get an error:

  File "/home/tuomas/sources/firedrake/src/firedrake/firedrake/parloops.py", line 210, in par_loop
    _map = _maps[measure.integral_type()]
KeyError: 'exterior_facet_top'

I therefore implemented the top/bottom loop with op2.par_loop.

[tkarna]

OK I think I got the top/bottom surfaces fixed using bt_masks. The lateral boundaries are still wrong.

[wence-]

This requires an extension in firedrake, see firedrakeproject/firedrake#751

[wence-]

Typo in test file name

[wence-]

UFL has symbolic operators for FacetAvg and CellAvg. We've never had a use case for them, so have not implemented them. However, we plausibly could.

I think we have a map from exterior facets to the nodes that are (topologically or geometrically) on the facet. However, it's only really used internally and I think it doesn't fit the bill.

You can do it like this, I think:

boundary_dofs = facet_support_dofs(V.fiat_element) # use vert_facet_support_dofs in extruded
local_facet_nodes = numpy.array(
            [boundary_dofs[e] for e in sorted(boundary_dofs.keys())])

kernel = """
void my_kernel(double **qmax, double **qmin, double **field, unsigned int *facet)
{
    local_facet_idx = 2d_array_of_local_facet_nodes;

    for (i = 0; i < nnode; i++ ) {
        idx = local_facet_idx[facet[0]][i];
        qmax[idx][0] = fmax(qmax[idx][0], field[idx][0]);
        ...
    }
}"""

op2.par_loop(my_kernel, V.mesh().exterior_facets.set, qmax.dat(op2.RW, qmax.exterior_facet_node_map()), qmin.dat(...), field.dat(...), V.mesh().exterior_facets.local_facet_dat(op2.READ))

Untested! This assumes all the fields are in the same space. If you have different spaces you'll need to encode the local_facet_idx array for each of them.

[tkarna]

Thanks. I got this far, just trying to write to the max field using the indexing. But it segfaults.

from FIAT.finite_element import facet_support_dofs
from FIAT.tensor_product import horiz_facet_support_dofs, vert_facet_support_dofs
import numpy

if self.is_2d:
    boundary_dofs = facet_support_dofs(self.P1CG.fiat_element)
    local_facet_nodes = numpy.array([boundary_dofs[e] for e in sorted(boundary_dofs.keys())])
    local_facet_idx = op2.Global(local_facet_nodes.shape, local_facet_nodes, dtype=np.int32, name='node_idx')

    code = """
    void my_kernel(double **qmax, double **qmin, double **field, unsigned int *facet, int **local_facet_idx)
    {
        for (int i = 0; i < %(nnodes)d; i++) {
            int idx = local_facet_idx[facet[0]][i];
            qmax[idx][0] = 4.5;
        }
    }"""
    facet_kernel = op2.Kernel(code % {'nnodes': local_facet_nodes.shape[0]}, 'my_kernel')

    op2.par_loop(facet_kernel,
                    self.P1CG.mesh().exterior_facets.set,
                    self.max_field.dat(op2.RW, self.max_field.exterior_facet_node_map()),
                    self.min_field.dat(op2.RW, self.min_field.exterior_facet_node_map()),
                    field.dat(op2.READ, field.exterior_facet_node_map()),
                    self.P1CG.mesh().exterior_facets.local_facet_dat(op2.READ),
                    local_facet_idx(op2.READ)
                    )

[tkarna]

For the time being I've omitted the lateral boundary treatment. I'd like to merge this version as the turbulence branch depends on it. We can revisit this issue this once migrated to new equation api.