ASPHER — Accelerated SPectral and HiERarchical contact solver

(pronounced "asper", as in asperity — or as the Latin asper, "rough")

ad astra per ASPHERa - through roughness, to the stars.

• Agentic framework: Frontier models (Anthropic & OpenAI)
• Coordinator: Vladislav A. Yastrebov (CNRS, Mines Paris - PSL)

C++17 boundary-integral solver for normal and frictional (tangential) contact of an elastic half-space — free-space or doubly periodic: exact Love/Boussinesq and Cerruti influence coefficients applied through hierarchical and spectral operators (H-matrix ACA, a matrix-free H²/bbFMM, exact FFT convolution), a |q| spectral preconditioner, nested-grid continuation and a mixed-precision Polonsky–Keer CG with certified stopping, with a Python interface.

The design combines an exact kernel with fast, matrix-free operators:

• Kernel: Love (1929) closed-form integration of the Boussinesq point solution over square elements — exact for every influence coefficient, served from an O(N) lookup table (translation invariance).
• Operator (default): matrix-free H²/FMM — black-box fast multipole method with Chebyshev interpolation (Fong & Darve 2009), all transfer and coupling operators cached by translation invariance.
• O(N) memory and matvec: 5.3 MiB at 512², ~13 B/DOF asymptotically; grids up to 16384² (2.7×10⁸ DOFs) run on a 32 GiB workstation.
• An exact FFT-convolution operator (backend="fft", zero-padded Love-kernel convolution) equals the dense matvec to roundoff (~10⁻¹⁵ rel L2) and is modestly faster than H² at Ns ≤ 2048 (~1.5×), ≈parity at Ns=4096; H² remains preferred for very large Ns.
• A doubly periodic operator (backend="periodic", spectral compliance 2/(E*|q|) on the Ns×Ns torus) solves periodic problems, validated against the exact Westergaard solution. A classical H-matrix (ACA) backend and a dense backend are kept for validation and small grids.
• Solver: Polonsky & Keer (1999) projected CG with overlap correction, accelerated by a |q| spectral preconditioner (a 13-tap real-space stencil on the contact set by default, or full-grid FFT), nested-grid (cascadic/FMG) continuation, an active-set mode that restricts the work to candidate contact points, warm starts, and an optional single-precision mode with double-accumulated reductions — up to ~15× faster than the cold double solve at equal contact-area accuracy. Convergence is certified on the returned pressure (Frank–Wolfe gap, penetration and load), and a solve that stagnates is reported as such, never as converged.
• Friction: Coulomb, Tresca or user-defined laws via Cerruti kernels, a projected CG for the tangential problem with force or displacement control, and a load-history driver (FrictionSolver).

Install

pip install aspher            # from PyPI (builds from the sdist)
# or from a checkout:
pip install .

Building from source needs a C++17 compiler and CMake ≥ 3.18; Eigen is found on the system or fetched automatically, and the FFT engine (pocketfft) is bundled. If the default compiler misbehaves (see the conda note in CLAUDE.md), override it with CMAKE_ARGS="-DCMAKE_CXX_COMPILER=/usr/bin/g++" pip install .. Wheels are built without -march=native (portable); local dev builds keep it.

Development build (see CLAUDE.md for the environment specifics):

cmake -S . -B build -DCMAKE_BUILD_TYPE=Release \
      -Dpybind11_DIR=$(python -m pybind11 --cmakedir)
cmake --build build -j && ctest --test-dir build

The Python module aspher*.so is placed in python/.

Release to PyPI (maintainers)

Bump version in pyproject.toml and CITATION.cff and commit, then build the sdist from a clean export of the commit — scikit-build-core packs every file in the directory that git does not ignore, so building in a working tree also ships untracked notes, local build directories and data:

python -m pip install -U build twine
rm -rf /tmp/aspher-release && mkdir /tmp/aspher-release
git archive HEAD | tar x -C /tmp/aspher-release
(cd /tmp && python -m build --sdist -o "$OLDPWD/dist" /tmp/aspher-release)
python -m twine check dist/aspher-X.Y.Z.tar.gz
python -m twine upload dist/aspher-X.Y.Z.tar.gz  # sdist only: PyPI rejects
                                                 # non-manylinux binary wheels

(python -m build is run from outside the checkout because the local build/ directory would shadow the build package.) The README is the PyPI project page, so review it before each release. The sdist is all-BSD (bundled pocketfft); binary wheels would need cibuildwheel (manylinux + auditwheel) and are a later step.

Quick start

import numpy as np
import rfgen          # random rough surfaces (installed with aspher)
import aspher as hc

def surface(N, rms=0.002):   # self-affine rough surface, Hurst 0.8
    h = rfgen.selfaffine_field(N=N, Hurst=0.8, k_low=4 / N, k_high=64 / N)
    return h * rms / h.std()

# a rigid rough indenter pressed onto an elastic half-space: gap0 = -height
Ns = 1024                              # power of two for the H2 operator
gap = -surface(Ns).ravel()

# one-call nested-grid solve (coarse -> fine, preconditioned, warm-started)
res = hc.solve_nested(grid_size=Ns, gap=gap, p_nominal=0.005,
                      coarsest=64, q=6)
print(res.status, res.contact_area, res.iterations)

# large grids: active set + single precision + light result
res = hc.solve_nested(grid_size=4096, gap=-surface(4096).ravel(),
                      p_nominal=0.005, active_set=True,
                      single_precision=True, light_result=True)

# or the explicit operator interface
# backend="periodic" is doubly periodic; h2/fft/hmatrix/dense are free-space
solver = hc.ContactSolver(grid_size=Ns, backend="h2", q=6)
res = solver.solve(gap, p_nominal=0.005, tol=1e-8, precond="fourier")
res.pressure        # (Ns, Ns), mean == p_nominal
res.contact_area    # Ac/A
res.fw_error, res.penetration_error   # the certified quantities (<= tol)

# doubly periodic contact (same API, also in solve_nested(backend="periodic"))
per = hc.ContactSolver(grid_size=Ns, backend="periodic")
res = per.solve(gap, p_nominal=0.005, tol=1e-8, precond="fourier-fft")

Frictional contact:

fs = hc.FrictionSolver(grid_size=256, E_star=1.0, nu=0.3,
                       model=hc.CoulombFriction(mu=0.3))
fs.set_gap(-surface(256))
fs.step(p_bar=0.05)                    # normal load
r = fs.step(q_bar=(0.01, 0.0), dt=1.0) # tangential force (total)
print(r.n_stick, r.n_slip, r.dissipation)

python example_rough_contact.py runs an end-to-end rough-surface demo (self-affine surface → applied pressure → contact map); python example_friction.py a Cattaneo–Mindlin partial-slip demo.

Performance highlights

Measured on a 20-core workstation (fixed-band self-affine roughness):

Quantity	Value
H² operator memory at Ns=2048 (4.2M DOFs)	51 MiB (dense: 128 TiB)
H² build / matvec at Ns=512	0.03 s / 9 ms (H-matrix: 24 s / 144 ms)
Nested solve Ns=1024, double, tol 1e-8	1.2 s
Nested solve Ns=4096, double, tol 1e-8	~40 s (50 it)
Nested solve Ns=4096, float + light result	~11 s (18 it, area matches double to 0.2%)
Ns=16384 (2.7×10⁸ DOFs), float, 32 GiB node	~25 GiB peak (standard path)
Ns=16384, active-set, float / double	125 s, 9.1 GiB / 315 s, 9.3 GiB (same contact area)
Periodic nested solve, Ns=8192, 8% contact, tol 1e-11	62 s on 10 threads / 161 s on 1 (49 it)
FFT matvec vs dense (rel L2, double / float)	~1×10⁻¹⁵ / ~1.4×10⁻⁷ (exact operator)
Matvec fft vs H² (q=6), Ns=1024 / 2048 / 4096	17 / 73 / 337 ms vs 27 / 110 / 331 ms (1.60× / 1.50× / 0.98×)*
Nested solve Ns=4096, p̄=0.002, fft vs h2, double	251 s (110 it) vs 177 s (91 it), areas agree to ~10⁻⁶*
Nested solve Ns=4096, p̄=0.002, fft vs h2, float	62 s (53 it) vs 51 s (52 it)*

* bench_fft.py, measured under desktop co-tenancy; ratios more reliable than absolutes. The fft/h2 nested iteration counts differ (both are valid PCG paths to the same solution).

Tests and benchmarks

• ctest --test-dir build — kernel vs analytics, H/H² matvec vs dense, Hertz and Westergaard contact vs theory, preconditioner symmetry and mixed-precision consistency, friction benchmarks (Cattaneo–Mindlin, Ciavarella–Jäger).
• python -m pytest tests/*_py.py tests/test_nested.py tests/test_h2.py — the Python API; both suites run in CI on Python 3.10–3.12.
• python bench_h2.py, bench_h2_memory.py, bench_h2_cputime.py — H² vs H-matrix and O(N) scaling sweeps up to Ns=16384.
• python bench_fft.py — FFT-convolution backend vs H²: matvec sweep (min-of-30) plus the four Ns=4096 nested solves (fft/h2 × double/float).
• python experiments/bench_cpp_precond.py — preconditioner + nested-grid acceleration (e.g. 180→45 iterations at Ns=1024).

License

The ASPHER source is BSD-3-Clause (see LICENSE).

The FFT transforms — the spectral preconditioner and the backend="fft" operator — run, by default, on the bundled pocketfft (BSD-3-Clause) — so the default build, including binary wheels, is permissively licensed throughout. Configuring with -DASPHER_USE_FFTW=ON switches to FFTW3 plans, which are slightly faster end-to-end (measured at Ns=4096: ~16% on a double solve, ~5% on a float solve, identical results) — but FFTW is GPL, so binaries built with that option are governed by the GPL as a whole. Choose per deployment:

CMAKE_ARGS="-DASPHER_USE_FFTW=ON" pip install .   # faster, GPL binaries
pip install .                                      # default, all-BSD

References

The methods implemented here are described in:

• J. Boussinesq, Application des potentiels à l'étude de l'équilibre et du mouvement des solides élastiques, Gauthier-Villars, 1885 — half-space point-load solution.
• A.E.H. Love, "The stress produced in a semi-infinite solid by pressure on part of the boundary", Phil. Trans. R. Soc. A 228 (1929) — exact element integrals.
• I.A. Polonsky, L.M. Keer, "A numerical method for solving rough contact problems based on the multi-level multi-summation and conjugate gradient techniques", Wear 231 (1999) 206–219 — the projected CG with overlap correction.
• R. Pohrt, Q. Li, "Complete boundary element formulation for normal and tangential contact problems", Phys. Mesomech. 17(4) (2014) 334–340 — the Cerruti element integrals for frictional contact.
• G. de Saxcé, Z.-Q. Feng, "The bipotential method: a constructive approach to design the complete contact law with friction and improved numerical algorithms", Math. Comput. Model. 28(4–8) (1998) 225–245 — the reference solver for frictional contact.
• W. Fong, E. Darve, "The black-box fast multipole method", J. Comput. Phys. 228 (2009) 8712–8725 — the Chebyshev H²/FMM operator.
• M. Bebendorf, Hierarchical Matrices, Springer, 2008; W. Hackbusch, Hierarchical Matrices: Algorithms and Analysis, Springer, 2015 — the H-matrix/ACA backend.
• V.A. Yastrebov, G. Anciaux, J.-F. Molinari, "From infinitesimal to full contact between rough surfaces: evolution of the contact area", Int. J. Solids Struct. 52 (2015) 83–102 — rough-contact context and validation methodology.

Solver theory (projected CG, spectral preconditioning, finite-precision implementation) is documented in doc/theory/pcg.tex; the H²/FMM operator in doc/theory/h2_fmm_detailed.tex; frictional contact in doc/theory/friction.tex.

Citing

If ASPHER contributes to your research, please cite it (see also CITATION.cff, which GitHub renders as a "Cite this repository" button):

@software{aspher,
  author  = {Yastrebov, Vladislav A. and {Claude (Anthropic)}},
  title   = {{ASPHER}: Accelerated {SP}ectral and {H}i{ER}archical contact solver},
  url     = {https://github.com/cofebem/aspher},
  version = {0.1.3.post1},
  year    = {2026},
}

Layout

include/, src/       kernels (Love, Cerruti), cluster tree, H-matrix, H2/FMM,
                     FFT and periodic operators, preconditioners, nested and
                     active-set solves, contact and friction solvers
python/bindings.cpp  pybind11 module `aspher`
third_party/         bundled pocketfft (BSD-3-Clause)
tests/               C++ tests (CTest) and Python tests (pytest)
.github/workflows/   CI: build + both test suites on Python 3.10-3.12
doc/theory/          PCG, H2/FMM and friction theory (LaTeX)