Complex `scale` implicitly converted to real

[conpierce8]

Motivation

The propagation of waves in periodic media is governed by the eigenvalue problem $$L u = \omega^2 u$$ where $L$ is a linear operator that is periodic with a period defined by the unit cell of the periodic medium. The eigenfunctions satisfy the Bloch form: $$u(x) = \hat{u}(x) e^{-i k \cdot x}$$ where $k$ is the wavevector and $\hat{u}(x)$ is periodic on the unit cell (i.e. $\hat{u}(x) = \hat{u}(x + a)$, where $a$ is a basis vector of the unit cell). Note that the eigenfunction $u(x)$ is Floquet-periodic on the unit cell, i.e. the solution on one boundary is related to the solution on the opposite boundary by $u(x) = u(x + a) e^{i k \cdot a}$. A common solution technique uses the finite element method to seek solutions $u(x)$ subject to Floquet periodic boundary conditions. For example, one might consider the unit square subject to the boundary conditions $$u(x=1) = u(x=0) e^{-i k_1}$$ where $k_1$ is the x-component of the wavenumber. This can be represented as a multi-point constraint with slave dofs located at $x=1$, master dofs located at $x=0$, and a (complex-valued) scale of $e^{-i k_1}$.

Issue

It appears that complex-valued scales are not fully implemented in dolfinx_mpc. Consider the following MWE, heavily adapted from demo_periodic_gep.py, which attempts to apply a periodic constraint with a complex-valued scale. A warning is issued suggesting that an implicit conversion from complex to real takes place, i.e. that the imaginary part of the scale is ignored. It appears that this error arises at the cpp wrapper layer, as create_periodic_constraint_topological takes a double scale, rather than a PetscScalar scale.

WARNING:py.warnings:/usr/local/lib/python3.10/dist-packages/dolfinx_mpc/multipointconstraint.py:110: ComplexWarning: Casting complex values to real discards the imaginary part
  mpc_data = dolfinx_mpc.cpp.mpc.create_periodic_constraint_topological(

It looks like there is an intention to support this, but perhaps it was never finished? mpc_data stores coefficients using the PETSc scalar type, but all of the call signatures in PeriodicConstraint accept only reals.

MWE

import dolfinx.fem as fem
import numpy as np
from dolfinx.mesh import create_unit_square, locate_entities_boundary, meshtags
from dolfinx_mpc import MultiPointConstraint
from mpi4py import MPI

N = 1
mesh = create_unit_square(MPI.COMM_WORLD, N, N)
V = fem.FunctionSpace(mesh, ("Lagrange", 1))
fdim = mesh.topology.dim - 1

def pbc_slave_to_master_map(x):
    out_x = x.copy()
    out_x[0] = x[0] - 1
    return out_x

def pbc_is_slave(x):
    return np.isclose(x[0], 1)

def pbc_is_master(x):
    return np.isclose(x[0], 0)

# Parse boundary conditions
bcs = []

pbc_slave_tag = 2
facets = locate_entities_boundary(mesh, fdim, pbc_is_slave)
arg_sort = np.argsort(facets)
pbc_meshtags = meshtags(mesh,
                        fdim,
                        facets[arg_sort],
                        np.full(len(facets), pbc_slave_tag, dtype=np.int32))
pbc_scale = np.exp(1j * np.pi / 5) # arbitrary complex-valued scale

# Create MultiPointConstraint object
mpc = MultiPointConstraint(V)
mpc.create_periodic_constraint_topological(V, pbc_meshtags,
                                           pbc_slave_tag,
                                           pbc_slave_to_master_map,
                                           bcs,
                                           pbc_scale)

Environment

The above was run in a complex-only Docker container that I built from the following recipe:

FROM dolfinx/dev-env:nightly as dolfinx-mpc

ARG PETSC_SCALAR_TYPE=complex
ARG DOLFINX_CMAKE_BUILD_TYPE=RelWithDebInfo

WORKDIR /tmp

# Set env variables
ENV PETSC_ARCH="linux-gnu-${PETSC_SCALAR_TYPE}-32" \
    HDF5_MPI="ON" \
    CC=mpicc \
    HDF5_DIR="/usr/local"

# Install DOLFINx dependencies:
#   - basix
#   - ufl
#   - ffcx
ADD basix /tmp/basix
ADD ufl /tmp/ufl
ADD ffcx /tmp/ffcx
RUN PYBIND11_DIR="$(python3 -c 'import sys, pybind11; sys.stdout.write(pybind11.get_cmake_dir())')" && \
    # pybind11 needs to be manually located first
    cd basix && \
    cmake -G Ninja -DCMAKE_BUILD_TYPE=${DOLFINX_CMAKE_BUILD_TYPE} "-Dpybind11_DIR=${PYBIND11_DIR}" -B build-dir -S .  && \
    cmake --build build-dir && \
    cmake --install build-dir && \
    pip3 install -v -e ./python  && \
    cd ../ufl && \
    pip3 install --no-cache-dir . && \
    cd ../ffcx && \
    pip3 install --no-cache-dir .

# Install dolfinx
ADD dolfinx /tmp/dolfinx
RUN cd dolfinx && \
    mkdir -p build && \
    cd build && \
    cmake -G Ninja -DCMAKE_BUILD_TYPE=${DOLFINX_CMAKE_BUILD_TYPE} ../cpp && \
    ninja install && \
    . /usr/local/lib/dolfinx/dolfinx.conf && \
    cd ../python && \
    pip3 install --no-dependencies --no-cache-dir .

# Install dolfinx_mpc
ADD dolfinx_mpc /tmp/dolfinx_mpc
RUN cd dolfinx_mpc && \
    mkdir -p build && \
    cd build && \
    cmake -G Ninja -DCMAKE_BUILD_TYPE="${DOLFINX_CMAKE_BUILD_TYPE}" ../cpp && \
    ninja install -j4 && \
    cd .. && \
    pip3 install -v --no-deps --upgrade python/.

# Install h5py
#RUN pip3 install --no-binary=h5py h5py meshio

WORKDIR /root

Library versions:

• basix: 1dd043332314a1633f3bdc32bb15ce31e498d892
• UFL: fd773880195b39327f652a47a892d89f7ef200cc
• ffcx: c007abd0753fb7a5a1f6c14ef6ff81b099666043
• dolfinx: 66f5adfc2d4e4c32e1e50a7ac472337b591ccb28
• dolfinx_mpc: 1a1520862e5ef2f05c0790a2ea94d2b231d01ff1

[conpierce8]

Closing this issue as complex values now appear to be treated correctly.