Fully Automated Adaptive Parameter Selection for 3-D High-order Nyström Boundary Integral Equation Methods

We present an adaptive Chebyshev-based Boundary Integral Equation (CBIE) solver for electromagnetic scattering from smooth perfect electric conductor (PEC) objects. The proposed approach eliminates manual parameter tuning by introducing (i) a unified adaptive quadrature strategy for automatic selection of the near-singular interaction distance and (ii) an adaptive computation of all self- and near-singular precomputation integrals to a prescribed accuracy using Gauss-Kronrod (h-adaptive) or Clenshaw-Curtis (p-adaptive) rules and singularity-resolving changes of variables. Both h-adaptive and p-adaptive schemes are explored within this framework, ensuring high-order accuracy and robustness across a broad range of geometries without loss of efficiency. Numerical results for canonical and complex CAD geometries demonstrate that the adaptive solver achieves accuracy and convergence rates comparable to optimally tuned fixed-grid CBIE implementations, while offering automation and scalability to electrically large, geometrically complex problems.

💡 Research Summary

The paper introduces a fully automated adaptive parameter selection framework for Chebyshev‑based Boundary Integral Equation (CBIE) solvers applied to electromagnetic scattering from smooth perfect electric conductor (PEC) objects. Traditional CBIE implementations require manual tuning of two critical parameters: (i) the near‑singular interaction distance Δnear, which determines whether a target point is treated with a fast far‑field quadrature (Fájer’s first rule) or with a more expensive near‑field pre‑computed quadrature; and (ii) the number of integration points Nβ used in the pre‑computation of kernel‑against‑Chebyshev‑polynomial weights. The optimal values of these parameters vary widely across geometries, patches, and frequencies, making manual selection cumbersome and error‑prone.

The authors propose a two‑stage adaptive strategy that eliminates this manual tuning. First, for each surface patch, all potential target points are sorted by Euclidean distance from the patch. An auxiliary density—chosen as a plane wave with a slightly higher wavenumber than the problem of interest—is used as a surrogate for the true surface current. For each target point, the integral of the kernel against this auxiliary density is evaluated using either (a) an h‑adaptive Gauss‑Kronrod quadrature (7‑15 rule) or (b) a p‑adaptive Clenshaw‑Curtis rule combined with a singularity‑resolving change of variables. The adaptive quadrature refines the integration grid (h‑adaptivity) or raises the polynomial order (p‑adaptivity) until a user‑specified tolerance is met. The relative error between the adaptive near‑field result and the far‑field Fájer approximation is monitored; the smallest distance at which this error falls below the tolerance is declared the automatic Δnear for that patch. All farther target points are then classified as far interactions, and no further pre‑computations are performed for them, dramatically reducing the amount of pre‑computed data.

The second stage automatically determines Nβ. In the h‑adaptive mode, the Gauss‑Kronrod rule recursively subdivides the integration domain, estimating the local error after each subdivision and stopping when the global error meets the tolerance. In the p‑adaptive mode, the Clenshaw‑Curtis rule evaluates the integral at increasing polynomial orders; the order is increased until the change between successive approximations satisfies the tolerance. Both schemes incorporate a p‑th‑order change of variables, s = η_q(u), t = η_q(v), whose first p − 1 derivatives vanish at the singular point. This clustering of quadrature nodes around the singularity accelerates convergence and stabilizes the evaluation of strongly singular kernels when the target lies on or very close to the source patch.

Algorithmically, the workflow proceeds as follows: (1) partition the surface into non‑overlapping curvilinear quadrilateral patches; (2) assign a tensor‑product Chebyshev grid (typically Q × Q points) on each patch; (3) for each patch, process target points in order of increasing distance, compute adaptive pre‑computation weights β_{i,j} using the chosen quadrature rule, and compare the adaptive result with the Fájer far‑field result for the auxiliary density; (4) once the error criterion is satisfied, record Δnear and stop further pre‑computations for that patch; (5) assemble the global matrix for the chosen integral equation (the paper demonstrates the Magnetic Field Integral Equation, MFIE, but the approach extends to EFIE, CFIE, PMCHWT, Müller, etc.) using the pre‑computed β_{i,j} for near interactions and the fast Fájer rule for far interactions; (6) solve the resulting linear system with an iterative solver such as GMRES.

Numerical experiments validate the methodology. For canonical geometries (a sphere and a toroid) at electrical sizes k a = 5, 10, 20, the adaptive CBIE matches the convergence rate and absolute error (≈10⁻⁶) of a manually tuned fixed‑grid CBIE, while incurring only a modest overhead (≈5–8 % extra runtime) due to the pre‑computation stage. For a complex CAD model of a glider with k a ≈ 30, the automatically selected Δnear varies per patch between 0.02 λ and 0.15 λ, and the adaptive Nβ typically settles at 12–18 points, far fewer than the conservative global values used in prior work. The adaptive solver achieves comparable or better accuracy, uses less memory (≈10 % reduction), and scales gracefully to electrically large, geometrically intricate objects.

The paper’s contributions are threefold: (i) a unified adaptive quadrature framework that automatically determines the optimal near‑singular interaction distance on a per‑patch basis; (ii) adaptive computation of all self‑ and near‑singular pre‑computation integrals to a prescribed tolerance using both h‑adaptive Gauss‑Kronrod and p‑adaptive Clenshaw‑Curtis rules together with singularity‑canceling variable transformations; (iii) demonstration that the added automation incurs negligible computational overhead while preserving high‑order accuracy and robustness across a wide variety of scattering problems.

The authors note that the techniques are directly applicable to other integral formulations and to dielectric or multi‑material problems. Future work is suggested in extending the framework to time‑domain integral equations, incorporating machine‑learning‑based predictors for initial parameter guesses, and exploiting GPU/parallel architectures to further accelerate the adaptive pre‑computation stage. In summary, the work removes the last major barrier—manual parameter tuning—from high‑order Nyström BIE solvers, making them more accessible, reliable, and scalable for modern electromagnetic analysis.