RFC: BlochReducer — a Bloch-Floquet periodic boundary-condition operator + phononic crystals examples (band structure, transmission)

[NicolasCayuela]

Hi,

I use COMSOL (via mph) for Python scripting in my research, and TensorMesh is
going to save me a lot of time — I'll be sure to cite the paper.

Following work on 2D phononic-crystal band-structure and transmission examples, I
would like to discuss whether TensorMesh should gain a small public operator for
Bloch-Floquet periodic boundary conditions for reasearch on metamaterials.

Motivation

Band structures and other Bloch-periodic eigenproblems (phononic/photonic
crystals, lattice dynamics, periodic acoustics) all need the same primitive:
reduce an assembled operator to the independent DOFs of a unit cell with a
wavevector-dependent Floquet phase,

A_r(k) = T(k)ᴴ A T(k), u(r + R) = exp(i k·R) u(r).

TensorMesh has no primitive for this today, so every periodic study hand-rolls
the node pairing and the phase reduction (I did this across the phononic examples
before factoring it out). A small operator removes that duplication and is the
natural periodic counterpart of Condenser (which eliminates Dirichlet DOFs).

Proposed scope

A single operator module in the existing operator family:

tensormesh/operator/bloch.py
    BlochReducer

exported as tensormesh.BlochReducer, documented in
docs/source/api/operator.rst next to Condenser, plus tests and examples
(below). I would not introduce a top-level tensormesh.bloch package.

Proposed API

from tensormesh import BlochReducer

bloch = BlochReducer(points, lattice_vectors, dofs_per_node=1)

Kr, Mr = bloch.reduce_system(K, M, k)   # T(k)^H K T(k), T(k)^H M T(k)
u_full = bloch.expand(u_reduced, k)     # scatter back with the Floquet phase

An nn.Module like Condenser (pairing buffers move with .to(device)).
Design:

• pairs opposite unit-cell faces by greedy lattice translation in float space, so
  non-orthogonal (e.g. triangular) lattices work;
• builds Tᴴ A T without forming T: remap A's COO indices onto the masters,
  weight by conj(φ_i) φ_j, coalesce — one pass over the non-zeros;
• scalar or vector DOFs (dofs_per_node); partial periodicity (1–3 lattice
  vectors, e.g. a waveguide);
• accepts a SparseMatrix (→ SparseMatrix) or a dense torch.Tensor (→ dense)
  for hand-assembled beam/truss systems;
• precondition: a periodic mesh (opposite faces carry matching nodes, e.g. gmsh
  setPeriodic) — same spirit as Condenser needing a boundary mask.

I am not attached to the exact name (BlochReducer / FloquetReducer / …) or to
the reduce/expand surface; I would prefer to match TensorMesh's existing
operator conventions.

Validation plan

• tests/operator/test_bloch.py (13 cases): node pairing (square, corner,
  non-orthogonal), reduction equals an explicit dense Tᴴ A T, Hermiticity of
  the reduced operator, k=0 (real, phase 1), partial periodicity, vector DOFs,
  dense/sparse equivalence;
• CPU + float64;
• clean docs build (autoclass in api/operator.rst);
• three worked examples validating the operator end-to-end against COMSOL
  Pressure-Acoustics: square band structure (rigid, ~0.1%), triangular
  band structure (penetrable two-medium, non-orthogonal, ~0.5%), finite-slab
  transmission (mean |ΔT| 0.007).

Examples

examples/wave/phononic_crystal_2d/
    band_structure_square.py
    band_structure_triangular.py
    transmission_slab.py

Example-only (no extra public API), depending only on
BlochReducer + the existing assemblers.

Square-lattice band structure — rigid cylinders in water, M-Γ-X-M (~0.1 % vs COMSOL):

Triangular-lattice band structure — penetrable steel in water, non-orthogonal lattice, M-Γ-K-M (~0.5 % vs COMSOL):

Transmission through a finite slab — 10 rigid cylinders, first-order radiation BC (mean |ΔT| 0.007 vs COMSOL):

Working code (operator + tests + examples + docs) is ready on a branch in my fork:
https://github.com/NicolasCayuela/TensorMesh/tree/bloch-floquet%2Bphononic_crystal

Feedback welcome

Mainly I'd like to know whether a Bloch-Floquet periodic operator is something
you'd want in TensorMesh proper. I'm happy to match your conventions — the name,
the file placement, and whether the surface should be reduce/expand or a
Condenser-style __call__ are all open. If it's of interest I'll send a PR.

Beyond this, I'd like to keep contributing — 3D examples and further work on
metamaterials for elastic and acoustic waves.

Thanks!

[Shizheng-Wen]

Hi Nicolas,

First, sorry for my late reply. I was overwhelmed by the teaching and research duties at the end of the semester.
Thanks for the your thorough RFC! and glad TensorMesh is useful for your metamaterials work!

Yes! This is something we'd want in TensorMesh! I'd be happy to take the PR and invite you as the core contributor for TensorMesh! A Bloch-Floquet reducer is a clean periodic counterpart to Condenser in the same operator/ family, and it lands right on our roadmap: complex-valued FEM was recently unblocked end-to-end (Helmholtz already validated), and acoustic/metamaterial problems are exactly the follow-up direction we named there. Your reduced operator is complex-Hermitian, so it also exercises the complex assembly + complex LDL^H solve path nicely.

Your proposed placement and scope match our conventions almost exactly, so please go ahead. A few small steers so the PR lands smoothly:

• Name & placement: BlochReducer in tensormesh/operator/bloch.py, exported as tensormesh.BlochReducer, autodoc'd next to Condenser in api/operator.rst — all good, keep the name.
• Surface: explicit reduce_system(K, M, k) / expand(u, k) is clearer than overloading __call__ here (k-dependent, two matrices), so I'd keep that. If you can, mirror Condenser's verb for the scatter-back (recover) or at least keep the nn.Module + buffers-move-with-.to() contract you already described.
• Vector DOFs: please document the component interleaving dofs_per_node assumes (node-major vs component-major) and make it match how the vector assemblers (e.g. linear elasticity) and the projector lay out DOFs — that's what will let the elastic-wave follow-ups drop in cleanly.
• Complex path: route the band-structure eig and any solves through torch-sla (eigsh for the generalized Hermitian problem, the complex solver for transmission) so the examples also cover the complex stack.
• Examples: examples/wave/phononic_crystal_2d/ is the right home (sibling to examples/wave/helmholtz/). Please make them runnable without COMSOL — commit the reference values as small arrays so CI/users can reproduce the ~0.1–0.5% comparison offline.
• Tests/CI: 13 cases on CPU+float64 is exactly our convention. The PR will need a DCO sign-off (git commit -s) on every
  commit and goes out under Apache-2.0 — see CONTRIBUTING.md. The merge gate is DCO + tests (py3.10/3.11/3.12).

One scoping suggestion: keep this first PR to the 2D operator + the three examples + tests + docs (You don't need to translate, just focus on the English version), and do 3D / elastic-wave extensions as follow-ups — easier to review, and we can add a gallery page once it's in. Looking forward to the PR, and great to have you contributing!

Best,
Shizheng

[NicolasCayuela]

Hi thanks for this kind message.

I created the PR. I'd like to continue contributing; I'm thinking of creating 3D examples in mechanics and acoustic resonators and maybe even some photonics examples as well when i have time.

Best,

Nicolas