empty-bdry.py

"""Reproduce bug that hits us in parallel."""

from mirgecom.mpi import mpi_entry_point
from functools import partial
import pyopencl as cl
import pyopencl.tools as cl_tools

from meshmode.array_context import PyOpenCLArrayContext
from meshmode.dof_array import thaw
from grudge.eager import EagerDGDiscretization
from grudge import sym as grudge_sym

from mirgecom.simutil import create_parallel_grid


def main(ctx_factory=cl.create_some_context):
    """Drive bug example."""
    cl_ctx = ctx_factory()
    queue = cl.CommandQueue(cl_ctx)
    actx = PyOpenCLArrayContext(queue,
                allocator=cl_tools.MemoryPool(cl_tools.ImmediateAllocator(queue)))

    dim = 3
    nel_1d = 17
    order = 1
    box_ll = -0.5
    box_ur = 0.5

    bnd_tags = {
        "left": ["-x", ],
        "right": ["+x", ],
        "bottom": ["-y", ],
        "top": ["+y", ],
        "back": ["-z", ],
        "front": ["+z", ],
    }

    boundaries = {
        grudge_sym.DTAG_BOUNDARY("left"): -1,
        grudge_sym.DTAG_BOUNDARY("right"): 1,
        grudge_sym.DTAG_BOUNDARY("bottom"): -2,
        grudge_sym.DTAG_BOUNDARY("top"): 2,
        grudge_sym.DTAG_BOUNDARY("back"): -3,
        grudge_sym.DTAG_BOUNDARY("front"): 3
    }

    num_parts = 4
    from meshmode.mesh.generation import generate_regular_rect_mesh
    mesh = generate_regular_rect_mesh(a=(box_ll,) * dim,
                            b=(box_ur,) * dim, n=(nel_1d,) * dim,
                            boundary_tag_to_face=bnd_tags)
    from meshmode.distributed import get_partition_by_pymetis
    part_per_element = get_partition_by_pymetis(mesh, num_parts)

    from meshmode.mesh.processing import partition_mesh
    local_mesh, _ = partition_mesh(mesh, part_per_element, 0)

    print(f"{local_mesh.nelements=},{mesh.nelements=}")

    from mpi4py import MPI
    discr = EagerDGDiscretization(actx, local_mesh, order=order)

    local_boundaries = {}
    nonlocal_boundaries = {}
    for btag in boundaries:
        bnd_discr = discr.discr_from_dd(btag)
        bnd_nodes = thaw(actx, bnd_discr.nodes())
        bnd_normal = thaw(actx, discr.normal(btag))
        num_bnd_nodes = bnd_nodes[0][0].shape[0]

        if num_bnd_nodes > 0:
            local_boundaries[btag] = boundaries[btag]
        else:
            nonlocal_boundaries[btag] = boundaries[btag]
            #print(f"{rank=},{btag=}")
            print(f"{bnd_nodes=}")
            print(f"{bnd_normal=}")

        # next line reproduces issue
        bnd_nodes @ bnd_normal
        #print(f"{rank=},{btag=},{result=}")

    #print(f"{rank=},{nonlocal_boundaries=}")


if __name__ == "__main__":
    main()