feat: add BEM radiation coupling via bempp-cl

main.nf

@@ -19,6 +19,9 @@ params.max_freq = 8000     // Maximum frequency for the sweep in Hz
 params.num_intervals = 100 // Number of frequency steps in the sweep
 params.mesh_size = 0.01    // Target mesh element size in meters
 
+// Radiation impedance model at the horn mouth
+params.radiation_model = "plane_wave"  // plane_wave, flanged_piston, unflanged_piston, bem
+
 // Execution Settings
 params.num_bands = 8       // Number of parallel jobs for the solver
 params.outdir = "./results"
@@ -79,7 +82,8 @@ process run_simulation {
         --max-freq ${max_f} \
         --num-intervals ${num_intervals_per_band} \
         --length ${params.length} \
-        --mesh-size ${params.mesh_size}
+        --mesh-size ${params.mesh_size} \
+        --radiation-model ${params.radiation_model}
     """
 }
 
@@ -227,7 +231,8 @@ process run_auto_simulation {
         --max-freq ${max_f} \
         --num-intervals ${num_intervals_per_band} \
         --length ${params.length} \
-        --mesh-size ${params.mesh_size}
+        --mesh-size ${params.mesh_size} \
+        --radiation-model ${params.radiation_model}
     """
 }
 

packages/horn-solver/Dockerfile

@@ -17,7 +17,8 @@ WORKDIR /app
 COPY ./packages /app/packages
 RUN pip install --no-deps /app/packages/horn-core \
     /app/packages/horn-drivers \
-    && pip install /app/packages/horn-solver
+    && pip install /app/packages/horn-solver \
+    && pip install bempp-cl
 
 # --- Test Stage ---
 FROM production as test

packages/horn-solver/scripts/compare_radiation_models.py

@@ -0,0 +1,191 @@
+#!/usr/bin/env python3
+"""Compare SPL curves from different radiation impedance models.
+
+Runs the horn solver with plane_wave, flanged_piston, and bem radiation
+models on the same STEP geometry, then plots an overlay comparison.
+
+Usage:
+    python compare_radiation_models.py \
+        --step-file horn.step \
+        --length 0.5 \
+        --min-freq 200 \
+        --max-freq 4000 \
+        --num-intervals 50 \
+        --output-dir comparison_results
+"""
+
+import argparse
+import sys
+from pathlib import Path
+
+import numpy as np
+import pandas as pd
+
+
+RADIATION_MODELS = ["plane_wave", "flanged_piston", "bem"]
+
+
+def run_model(
+    step_file: str,
+    length: float,
+    freq_range: tuple,
+    num_intervals: int,
+    mesh_size: float,
+    radiation_model: str,
+    output_dir: Path,
+) -> pd.DataFrame:
+    """Run the solver for a single radiation model and return results."""
+    from horn_solver.solver import run_simulation_from_step
+
+    output_file = output_dir / f"results_{radiation_model}.csv"
+    driver_params = {"length": length}
+
+    try:
+        run_simulation_from_step(
+            step_file=step_file,
+            freq_range=freq_range,
+            num_intervals=num_intervals,
+            driver_params=driver_params,
+            output_file=str(output_file),
+            max_freq_mesh=freq_range[1],
+            mesh_size=mesh_size,
+            radiation_model=radiation_model,
+        )
+        return pd.read_csv(output_file)
+    except (ImportError, RuntimeError) as exc:
+        print(f"WARNING: {radiation_model} failed: {exc}")
+        return None
+
+
+def plot_comparison(all_results: dict, output_dir: Path):
+    """Plot overlaid SPL curves for all models."""
+    try:
+        import matplotlib
+        matplotlib.use("Agg")
+        import matplotlib.pyplot as plt
+    except ImportError:
+        print("matplotlib not available — skipping plot")
+        return
+
+    fig, ax = plt.subplots(figsize=(10, 6))
+
+    colors = {"plane_wave": "blue", "flanged_piston": "orange", "bem": "red"}
+    labels = {
+        "plane_wave": "Plane Wave (Z=rho*c)",
+        "flanged_piston": "Flanged Piston",
+        "bem": "BEM (nonlocal)",
+    }
+
+    for model, df in all_results.items():
+        ax.semilogx(
+            df["frequency"], df["spl"],
+            color=colors.get(model, "gray"),
+            label=labels.get(model, model),
+            linewidth=1.5,
+        )
+
+    ax.set_xlabel("Frequency (Hz)")
+    ax.set_ylabel("SPL (dB)")
+    ax.set_title("Radiation Model Comparison")
+    ax.legend()
+    ax.grid(True, alpha=0.3)
+
+    plot_path = output_dir / "radiation_model_comparison.png"
+    fig.savefig(plot_path, dpi=150, bbox_inches="tight")
+    plt.close(fig)
+    print(f"Plot saved: {plot_path}")
+
+
+def compute_deviations(all_results: dict, output_dir: Path):
+    """Compute pairwise deviations between models and save to CSV."""
+    models = list(all_results.keys())
+    rows = []
+
+    for i in range(len(models)):
+        for j in range(i + 1, len(models)):
+            m1, m2 = models[i], models[j]
+            df1, df2 = all_results[m1], all_results[m2]
+
+            # Interpolate to common frequencies if needed
+            if np.allclose(df1["frequency"].values, df2["frequency"].values):
+                diff = np.abs(df1["spl"].values - df2["spl"].values)
+            else:
+                from scipy.interpolate import interp1d
+                f_common = np.union1d(df1["frequency"].values, df2["frequency"].values)
+                spl1 = interp1d(df1["frequency"], df1["spl"], fill_value="extrapolate")(f_common)
+                spl2 = interp1d(df2["frequency"], df2["spl"], fill_value="extrapolate")(f_common)
+                diff = np.abs(spl1 - spl2)
+
+            rows.append({
+                "model_1": m1,
+                "model_2": m2,
+                "max_deviation_dB": float(np.max(diff)),
+                "mean_deviation_dB": float(np.mean(diff)),
+                "std_deviation_dB": float(np.std(diff)),
+            })
+            print(f"  {m1} vs {m2}: max={np.max(diff):.2f} dB, mean={np.mean(diff):.2f} dB")
+
+    dev_df = pd.DataFrame(rows)
+    dev_path = output_dir / "model_deviations.csv"
+    dev_df.to_csv(dev_path, index=False)
+    print(f"Deviations saved: {dev_path}")
+
+
+def main():
+    parser = argparse.ArgumentParser(description="Compare radiation impedance models.")
+    parser.add_argument("--step-file", type=str, required=True)
+    parser.add_argument("--length", type=float, required=True)
+    parser.add_argument("--min-freq", type=float, default=200.0)
+    parser.add_argument("--max-freq", type=float, default=4000.0)
+    parser.add_argument("--num-intervals", type=int, default=50)
+    parser.add_argument("--mesh-size", type=float, default=0.01)
+    parser.add_argument("--output-dir", type=str, default="comparison_results")
+    parser.add_argument(
+        "--models", type=str, nargs="+", default=RADIATION_MODELS,
+        help="Radiation models to compare (default: all)",
+    )
+    args = parser.parse_args()
+
+    output_dir = Path(args.output_dir)
+    output_dir.mkdir(parents=True, exist_ok=True)
+
+    freq_range = (args.min_freq, args.max_freq)
+    all_results = {}
+
+    for model in args.models:
+        print(f"\n{'='*60}")
+        print(f"Running: {model}")
+        print(f"{'='*60}")
+
+        df = run_model(
+            step_file=args.step_file,
+            length=args.length,
+            freq_range=freq_range,
+            num_intervals=args.num_intervals,
+            mesh_size=args.mesh_size,
+            radiation_model=model,
+            output_dir=output_dir,
+        )
+        if df is not None:
+            all_results[model] = df
+
+    if len(all_results) < 2:
+        print("ERROR: Need at least 2 successful models to compare")
+        return 1
+
+    print(f"\n{'='*60}")
+    print("Computing deviations...")
+    print(f"{'='*60}")
+    compute_deviations(all_results, output_dir)
+
+    print(f"\n{'='*60}")
+    print("Generating comparison plot...")
+    print(f"{'='*60}")
+    plot_comparison(all_results, output_dir)
+
+    print(f"\nDone. Results in: {output_dir}")
+    return 0
+
+
+if __name__ == "__main__":
+    sys.exit(main())

packages/horn-solver/scripts/test_bempp_install.py

@@ -0,0 +1,63 @@
+#!/usr/bin/env python3
+"""Phase 1a: Verify bempp-cl is importable with the Numba backend.
+
+Run inside the dolfinx/dolfinx:v0.8.0 container after `pip install bempp-cl`:
+
+    docker run --rm dolfinx/dolfinx:v0.8.0 bash -c \
+        "pip install bempp-cl && python3 scripts/test_bempp_install.py"
+"""
+
+import sys
+
+
+def main() -> int:
+    # 1. Import bempp-cl (v0.4.x module name)
+    try:
+        import bempp.api as bempp_api
+    except ImportError:
+        print("FAIL: could not import bempp.api — is bempp-cl installed?")
+        return 1
+
+    print(f"OK: imported bempp.api (version {bempp_api.__version__})")
+
+    # 2. Verify Numba backend is active (no OpenCL)
+    device = getattr(bempp_api, "DEFAULT_DEVICE_INTERFACE", "unknown")
+    print(f"  DEFAULT_DEVICE_INTERFACE = {device!r}")
+    if device != "numba":
+        print(f"WARN: expected 'numba' backend, got {device!r}")
+        # Not fatal — the OpenCL backend works too, but Numba is the goal
+
+    # 3. Create a simple sphere grid and assemble a single-layer operator
+    try:
+        grid = bempp_api.shapes.regular_sphere(3)  # refinement level 3
+        print(f"OK: created sphere grid with {grid.number_of_elements} elements")
+    except Exception as exc:
+        print(f"FAIL: could not create sphere grid: {exc}")
+        return 1
+
+    try:
+        space = bempp_api.function_space(grid, "P", 1)
+        print(f"OK: created P1 function space ({space.global_dof_count} DOFs)")
+    except Exception as exc:
+        print(f"FAIL: could not create function space: {exc}")
+        return 1
+
+    # Assemble Helmholtz single-layer operator at k=1
+    k = 1.0
+    try:
+        slp = bempp_api.operators.boundary.helmholtz.single_layer(
+            space, space, space, k
+        )
+        mat = slp.weak_form()
+        print(f"OK: assembled Helmholtz single-layer operator (k={k})")
+        print(f"  Matrix shape: {mat.shape}")
+    except Exception as exc:
+        print(f"FAIL: could not assemble BEM operator: {exc}")
+        return 1
+
+    print("\nAll checks passed — bempp-cl is working with the Numba backend.")
+    return 0
+
+
+if __name__ == "__main__":
+    sys.exit(main())

packages/horn-solver/scripts/test_fem_bem_coupling.py

@@ -0,0 +1,97 @@
+#!/usr/bin/env python3
+"""Phase 1b: FEM-BEM coupling smoke test adapted for DOLFINx v0.8 + bempp-cl.
+
+Solves a simple exterior Helmholtz problem on a unit cube:
+  - FEM interior (DOLFINx) + BEM exterior (bempp-cl)
+  - Verifies the FEniCSx <-> bempp trace coupling machinery works
+
+Based on the mscroggs FEM-BEM coupling tutorial, updated for:
+  - DOLFINx v0.8 API (dolfinx.fem.functionspace, etc.)
+  - bempp-cl (bempp.api, not bempp_cl.api)
+
+Run inside the dolfinx container:
+    docker run --rm -v $PWD:/app dolfinx/dolfinx:v0.8.0 bash -c \
+        "source /usr/local/bin/dolfinx-complex-mode && \
+         pip install bempp-cl && python3 /app/packages/horn-solver/scripts/test_fem_bem_coupling.py"
+"""
+
+import sys
+import numpy as np
+
+
+def main() -> int:
+    # --- Imports ---
+    try:
+        import dolfinx
+        from dolfinx import fem, mesh
+        from mpi4py import MPI
+        import ufl
+        print(f"OK: DOLFINx {dolfinx.__version__}")
+    except ImportError as exc:
+        print(f"FAIL: DOLFINx import error: {exc}")
+        return 1
+
+    try:
+        import bempp.api as bempp_api
+        print(f"OK: bempp.api {bempp_api.__version__}")
+    except ImportError as exc:
+        print(f"FAIL: bempp import error: {exc}")
+        return 1
+
+    try:
+        from bempp.api.external import fenicsx as bempp_fenicsx
+        print("OK: imported bempp.api.external.fenicsx coupling module")
+    except ImportError as exc:
+        print(f"FAIL: bempp-fenicsx coupling not available: {exc}")
+        return 1
+
+    # --- Create a unit cube mesh in DOLFINx ---
+    print("\nCreating unit cube FEM mesh...")
+    domain = mesh.create_unit_cube(
+        MPI.COMM_WORLD, 5, 5, 5, cell_type=mesh.CellType.tetrahedron
+    )
+    V = fem.functionspace(domain, ("Lagrange", 1))
+    print(f"  FEM DOFs: {V.dofmap.index_map.size_global}")
+
+    # --- Extract boundary trace space ---
+    print("Extracting boundary trace space...")
+    try:
+        fenics_space, trace_matrix = bempp_fenicsx.fenics_to_bempp_trace_data(V)
+        print(f"  BEM trace DOFs: {fenics_space.global_dof_count}")
+        print(f"  Trace matrix shape: {trace_matrix.shape}")
+    except Exception as exc:
+        print(f"FAIL: trace extraction failed: {exc}")
+        return 1
+
+    # --- Assemble BEM operators ---
+    k = 1.0  # wavenumber
+    print(f"\nAssembling BEM operators (k={k})...")
+    try:
+        slp = bempp_api.operators.boundary.helmholtz.single_layer(
+            fenics_space, fenics_space, fenics_space, k
+        )
+        dlp = bempp_api.operators.boundary.helmholtz.double_layer(
+            fenics_space, fenics_space, fenics_space, k
+        )
+        print("  OK: single-layer and double-layer operators assembled")
+    except Exception as exc:
+        print(f"FAIL: BEM operator assembly failed: {exc}")
+        return 1
+
+    # --- Quick validation: apply operator to a constant function ---
+    print("Applying BEM operators to a test function...")
+    try:
+        ones = bempp_api.GridFunction(fenics_space, coefficients=np.ones(fenics_space.global_dof_count))
+        result = slp * ones
+        print(f"  OK: SLP * ones -> grid function with {result.coefficients.shape[0]} coefficients")
+        print(f"  Max coefficient magnitude: {np.max(np.abs(result.coefficients)):.6f}")
+    except Exception as exc:
+        print(f"FAIL: operator application failed: {exc}")
+        return 1
+
+    print("\nAll FEM-BEM coupling checks passed.")
+    return 0
+
+
+if __name__ == "__main__":
+    sys.exit(main())