Skip to content

Add additional nonlinear solvers and prepare to run polyfem with hybrid solver - #120

Draft
maxpaik16 wants to merge 24 commits into
polyfem:mainfrom
maxpaik16:max-dev
Draft

maxpaik16 wants to merge 24 commits into
polyfem:mainfrom
maxpaik16:max-dev

Conversation

@maxpaik16

Copy link
Copy Markdown
Contributor

No description provided.

maxpaik16 and others added 23 commits September 9, 2026 14:46
…em#114)

* Run the hybrid solver on hypre's threads-as-ranks backend

The hybrid solver needed MPI, which meant the whole application had to be
launched under mpirun. That is fine for a driver but not for a library:
anything embedding polysolve inherited the requirement.

hypre can now build its MPI surface on threads of one process
(HYPRE_ENABLE_THREAD_MPI), so the ranks become threads and mpirun goes
away. The resulting binary links no libmpi at all.

What this needed:

- tmpi_mpi_compat.hpp: the few MPI facilities hypre's backend does not
  already provide. hypre maps MPI_Bcast/Allreduce/Scatterv/... onto its
  own surface already; missing were the MPI-3 shared-memory windows
  (with threads, "shared memory" is just a pointer broadcast),
  MPI_IN_PLACE, and MPI_Init/Initialized/Finalized.

- The ranks are started by the solver rather than by mpirun. The first
  hybrid solver constructed calls hypre_tmpi_team_start(), the calling
  thread becomes rank 0 and keeps driving, and the workers sit in the
  existing command loop. They now return from their thread function
  instead of std::exit(0), and the last solver destroyed shuts the team
  down so the rest of the process sees a one-rank world again.

- is_running_worker_loop and worker_registry are thread_local. They were
  per-rank only because each rank used to be a process; as threads they
  were shared and the workers raced on them.

MPI remains available: POLYSOLVE_WITH_MPI still selects a rank-parallel
hybrid solver, it just no longer implies OpenMPI.

Co-Authored-By: Claude Opus 5 <noreply@anthropic.com>

* Add an SPMD benchmark driver and reporting scripts for the hybrid solver

The Catch2 tests cannot time the multi-rank path: they drive the solver
from rank 0 while the other ranks sit in the worker loop, and under
mpirun at more than one rank the suite is torn down when a worker leaves
the loop. Timing that measures the teardown, not the solve.

tests/bench_spmd.cpp is the same source built into both an OpenMPI and a
thread-MPI build, so the two backends are compared on one driver rather
than on two that happen to look similar. It builds a 7-point Laplacian
of a given grid size, solves it, and prints setup/solve/total plus the
relative residual so a run that converged differently cannot be mistaken
for a run that was merely faster.

scripts/bench_hybrid.sh sweeps rank counts on both backends. It takes an
flock so two sweeps cannot halve each other's scores, refuses to start
against a loaded machine rather than quietly reporting contended
numbers, and interleaves the backends inside each repetition so drift
hits both equally.

scripts/bench_report.py renders the CSV as Markdown: best-of-N rather
than mean, the worst run-to-run spread so the reader can judge whether a
gap is real, and a check that both backends reached the same residual at
each rank count.

Co-Authored-By: Claude Opus 5 <noreply@anthropic.com>

* Take nano-mpi as the MPI, instead of a backend inside hypre

The threads-as-ranks backend this branch was written against lived inside
hypre, behind a HYPRE_ENABLE_THREAD_MPI build mode, and needed a 183-line
compatibility shim here to fill the gaps. It is now a library of its own:

  https://github.com/danielepanozzo/nano-mpi

which changes what polysolve has to do. nano-mpi installs a header named
mpi.h, so hypre is built as an ordinary MPI build -- HYPRE_ENABLE_MPI=ON,
find_package(MPI) as it always did -- and the shim deletes entirely:

  * MPI_IN_PLACE, MPI_Init/Initialized/Finalized, and const-correct
    collectives are all native now.
  * MPI-3 shared-memory windows are native too. They were the last thing
    keeping the shim alive, and they are the one part of MPI's one-sided
    chapter that is trivially true when ranks are threads: the window is a
    real allocation plus everyone's offset into it. nano-mpi implements the
    layout MPI actually promises rather than the shim's simplification, so
    MPI_Win_shared_query answers correctly for every rank, not just rank 0.

  hypre_tmpi_team_start/join  ->  nanompi_team_start/join
  HYPRE_TMPI_NUM_THREADS      ->  NANOMPI_NUM_RANKS

The one piece of machinery this adds is cmake/nanompi-as-mpi/FindMPI.cmake.
hypre calls find_package(MPI REQUIRED) internally and would otherwise find a
system Open MPI -- whose ranks are processes, needing a launcher polysolve has
no way to invoke. The shim is on CMAKE_MODULE_PATH only when
POLYSOLVE_WITH_MPI is on, and it points MPI at nano-mpi.

It gives MPI::MPI_C the header path and nothing else, deliberately: hypre puts
that target in CMAKE_REQUIRED_LIBRARIES and runs check_c_source_compiles, and
try_compile() exports it into a scratch project where a link interface naming
nanompi::nanompi is a hard error. Whoever links MPI links nanompi::nanompi
explicitly instead. The probe in question is for MPI_Comm_f2c, which nano-mpi
does not declare -- it has no Fortran bindings and cannot, since Fortran SAVE
storage is per-process by language rule -- so it correctly comes out false.

Verified on macOS/arm64, Release:

  unit_tests: All tests passed (1521 assertions in 23 test cases)
  otool -L bench_spmd | grep -i mpi  ->  nothing

  bench_spmd 40 CPUHybrid, 64k unknowns:
    1 rank    setup 0.040  solve 0.098  relres 3.02e-11
    2 ranks   setup 0.028  solve 0.057  relres 7.57e-11
    4 ranks   setup 0.014  solve 0.041  relres 1.97e-11

nano-mpi is pinned to v0.1.0. hypre still points at the fork branch, because
the re-entrancy fix it needs (nine globals made thread-local) is still an open
PR upstream; that pointer moves when it lands.

Co-Authored-By: Claude Opus 5 <noreply@anthropic.com>

* Add a CI job that actually builds the hybrid solver

Every existing job configures with POLYSOLVE_WITH_MPI=OFF, which leaves the
hybrid solver out of the build entirely. A green CI on this branch therefore
said nothing about the thing the branch is for.

This job turns MPI on. Ranks are threads of the test process, so there is
nothing to install and nothing to launch -- the Linux dependency list is the
same as the others minus the mpi package.

Beyond building and running ctest it checks two things the port could plausibly
get wrong and still pass tests:

  * that no MPI runtime is linked into the binary, since the whole point is
    that ranks are threads and not processes;
  * that the solver converges at 1, 2 and 4 ranks, since a rank count that
    silently produces a wrong decomposition would otherwise look like a pass.

Co-Authored-By: Claude Opus 5 <noreply@anthropic.com>

* Move to nano-mpi v0.1.1

The Linux CI job added in the previous commit found what the macOS build did
not: AMGCL's MPI headers need MPI_CXX_DOUBLE_COMPLEX and MPI_CXX_FLOAT_COMPLEX
(amgcl/mpi/util.hpp maps std::complex<T> onto them unconditionally),
MPI_Exscan (mpi/partition/util.hpp) and MPI_Ialltoall
(mpi/coarsening/pmis.hpp). v0.1.1 has all four, plus the rest of the
nonblocking collectives and a Windows port.

Co-Authored-By: Claude Opus 5 <noreply@anthropic.com>

---------

Co-authored-by: Claude Opus 5 <noreply@anthropic.com>
NewtonCG is an inexact Newton method restricted to the CPU/GPU hybrid
linear solvers: it sets the linear solve's relative tolerance to
min(0.5, sqrt(||grad||)) each iteration and always accepts the
resulting (possibly inexact) direction.

NonlinearCG adds a first-order nonlinear conjugate gradient descent
strategy (Fletcher-Reeves, Polak-Ribiere, or Hestenes-Stiefel beta,
with negative-beta and periodic restarts to steepest descent).

Both are exposed as top-level solver choices and as entries in the
solver fallback list, with corresponding json spec options.

Co-Authored-By: Claude Sonnet 5 <noreply@anthropic.com>
AMGCL's recipe sets CMAKE_FIND_PACKAGE_PREFER_CONFIG before pulling AMGCL
in, so its own optional find_package(MPI) call skips Module mode (and
nano-mpi's FindMPI shim on CMAKE_MODULE_PATH) and resolves via Config mode
to a real system MPI instead. Since MPI::MPI_C is one global imported
target, CMake's bundled FindMPI unconditionally overwrites its
INTERFACE_INCLUDE_DIRECTORIES when that happens, silently re-pointing
hypre's #include <mpi.h> at a real MPI nano-mpi cannot run against and
crashing the CPUHybrid/NewtonCG solver with "reduction on datatype ... is
not supported" as soon as hypre assembles a real (multi-row) matrix.

Disable AMGCL's own MPI probe when nano-mpi is in play (its distributed
backend needs a real launcher nano-mpi doesn't provide anyway), and
promote CMAKE_MODULE_PATH to CACHE so the shim stays visible to directory
scopes that don't inherit the list(APPEND) through add_subdirectory().
AMGCL, CPUHybridSolver, GPUHybridSolver, HypreSolver, and MASSolver
previously read their block size (block_size/block_dim/dimension) out
of the JSON params inside set_parameters. Each now exposes a
set_block_size override instead, and set_parameters no longer reads
that key (removed from linear-solver-spec.json accordingly).
CPUHybridSolver's override broadcasts over MPI like set_parameters
does, so worker ranks stay in sync.

Solver::create(json, logger, strict_validation) gained a dimension
parameter that calls set_block_size before set_parameters. This is
threaded through polysolve::nonlinear (Problem::dimension(), Newton,
BFGS, and nonlinear::Solver::create) so a nonlinear solver's problem
dimension can reach its linear solver's block size.

Co-Authored-By: Claude Sonnet 5 <noreply@anthropic.com>
Wires set_block_mapping(VectorXi) through to HYPRE_BoomerAMGSetDofFunc
for HypreSolver, CPUHybridSolver, and GPUHybridSolver, so callers can
give BoomerAMG an explicit per-row function assignment instead of
relying on its default interleaved (row i -> i % dim) mapping.

HYPRE takes ownership of the dof_func array it's given, so it is
allocated with hypre_CTAlloc (matching the allocator HYPRE frees it
with) rather than new[]/malloc. GPUHybridSolver runs with
HYPRE_MEMORY_DEVICE, so the array is built on host and copied to a
HYPRE-owned device buffer via hypre_TMemcpy. CPUHybridSolver splits
rows across MPI ranks, and HYPRE_BoomerAMGSetDofFunc expects a
per-rank local array, so set_block_mapping broadcasts the full
mapping to every rank (mirroring the existing set_parameters/
set_block_size broadcast pattern) and each rank slices out its own
local rows once factorize() has partitioned them.

Adds a block_mapping test that checks an explicit mapping reproducing
HYPRE's own default interleaved assignment converges in the same
number of iterations as never calling set_block_mapping, for both
Hypre and CPUHybrid.

Co-Authored-By: Claude Sonnet 5 <noreply@anthropic.com>
Problem gains block_mapping(), returning an empty VectorXi by default,
that a problem subclass can override to report an explicit per-row
function (block) assignment. Newton (both sparse and dense solves)
and BFGS now check objFunc.block_mapping() before factorizing and
call linear_solver->set_block_mapping() when the problem provides
one, so a problem-specific mapping reaches the underlying multigrid
solver without callers having to wire it through manually.

Co-Authored-By: Claude Sonnet 5 <noreply@anthropic.com>
Lets a Problem report DOFs it considers problematic (e.g. poorly
conditioned or in contact) and contact patches, which Newton now
passes through to the linear solver each iteration. This is the
infrastructure AMGF and contact-patch Schwarz subdomain selection
build on in the hybrid solvers.

Co-Authored-By: Claude Sonnet 5 <noreply@anthropic.com>
Ports the following from the hybrid-dev branch onto the current
CPUHybridSolver/GPUHybridSolver (adapted to their nanompi worker-loop
and block-structure-based block-size handling, which hybrid-dev never
merged):

- select_bad_dofs_from_l1_row_norm: when disabled, the problematic
  subspace comes only from set_problematic_dofs (AMGF mode) instead of
  the row-norm heuristic, for both CPU and GPU Hybrid.
- subdomain_selection_strategy: KNEE, GMM (default, existing
  behavior), FD, and COST alternatives for thresholding row norms into
  the problematic subspace, for both CPU and GPU Hybrid.
  APOSTERIORI is reserved but rejected as not yet implemented.
- contact_patch_schwarz (CPU Hybrid only): gives each contact patch
  its own overlapping subdomain, accumulated additively (Schwarz) with
  the rest of the problematic subspace. Combined with
  select_bad_dofs_from_l1_row_norm=false this is AMGF-Schwarz.

Co-Authored-By: Claude Sonnet 5 <noreply@anthropic.com>
Adds select_bad_dofs_from_l1_norm, subdomain_selection_strategy, and
(CPU only) contact_patch_schwarz to the CPUHybrid/GPUHybrid spec
entries so the new options are discoverable and validated.

Co-Authored-By: Claude Sonnet 5 <noreply@anthropic.com>
Ported from the hybrid-dev branch's tests/linear_solve.cpp: load a
Matrix Market matrix (and optional rhs) and run any solver from a
JSON config, with warmup/repeat timing and a residual/peak-memory
report -- useful for reproducing or benchmarking a solve outside the
full polyfem integration.

Adapted to this branch's Solver API: dropped the -p/-e
positions/elements flags and the set_positions/set_elements calls
(Ichol-only, not part of this branch's Solver interface), and
dropped the MPI_Init/HYPRE_Initialize/rank-branching worker-loop
boilerplate, which was written for hybrid-dev's pre-nanompi,
real-MPI-process architecture -- CPUHybridSolver now initializes
HYPRE and spawns its own nano-mpi worker threads internally, so a
driver program needs none of that. -i/-t/--less_than (problematic
DOFs) are kept since set_problematic_dofs is part of this branch's
Solver API.

Co-Authored-By: Claude Sonnet 5 <noreply@anthropic.com>
NLOHMANN_JSON_SERIALIZE_ENUM(NormType, ...) was defined in Solver.cpp,
so its custom string<->enum mapping was only visible via ADL within
that translation unit. Any other caller of .get<NormType>() (e.g.
polyfem's NonlinearElasticVarForm, building the augmented-Lagrangian
sub-solver) silently fell back to nlohmann's generic enum-as-integer
conversion, which throws since json-specs declares norm_type as a
string ("L2", "Linf", "Euclidean"). Move the macro next to the enum
declaration in Problem.hpp instead; it expands to template<> functions,
so defining it in a header is safe.

Co-Authored-By: Claude Sonnet 5 <noreply@anthropic.com>
Log phase timers (analyze_pattern/factorize/solve, and the cuDSS/GPU
equivalents) at INFO instead of TRACE so scaling-experiment's
run_scaling_experiment.py can parse them out of default-level logs.
Route CPUHybridSolver's own phase logging through a CPUHYBRID_LOG_INFO
macro that only logs from rank 0 (previously rank!=0 workers were
silenced via spdlog::set_level(off), but that mutates the process-wide
default logger and also silences rank 0). Also reuse a single solver
instance across all -w/-r iterations in tests/linear_solve.cpp instead
of constructing a new one per iteration, since CPUHybrid's construction
spins up a multi-thread nanompi rank team and repeatedly tearing that
down and rebuilding it back-to-back invites a race in team create/teardown.

Co-Authored-By: Claude Sonnet 5 <noreply@anthropic.com>
CPUHybridSolver/GPUHybridSolver::name() reports "CPUAMGF"/"GPUAMGF"
instead of "CPUHybrid"/"GPUHybrid" whenever select_bad_dofs_from_l1_norm
is false (the informed/contact-patch subdomain-selection mode used by
AMGF and AMGF-Schwarz), but NewtonCG's constructor only accepted the
literal "CPUHybrid"/"GPUHybrid" names, rejecting the exact same solver
class under its AMGF display name.

Co-Authored-By: Claude Sonnet 5 <noreply@anthropic.com>
Callers that need a small, auxiliary linear solve unrelated to the
main Hessian solve (e.g. a constraint-elimination projection) have no
config-driven way to avoid inheriting a heavyweight solver (e.g.
Hybrid/AMGF) meant for the full-size system. Expose a simple_linear
string option, defaulting to Eigen::SimplicialLDLT, for such callers
to read instead of hardcoding the choice.
Both were left uninitialized as instance members and only ever set
inside pcg_solve()/the solve-info bookkeeping, so a solver queried via
get_info() before its first solve (or one whose pcg_solve exits before
reaching those assignments) could report garbage. Default them to 0.

Also drops a stray leftover empty statement.
@codecov

codecov Bot commented Sep 16, 2026

Copy link
Copy Markdown

Codecov Report

❌ Patch coverage is 51.13636% with 172 lines in your changes missing coverage. Please review.
✅ Project coverage is 75.95%. Comparing base (da4e7fe) to head (d6a9190).

Files with missing lines Patch % Lines
src/polysolve/linear/CPUHybridSolver.cpp 41.29% 118 Missing ⚠️
.../polysolve/nonlinear/descent_strategies/Newton.cpp 46.87% 17 Missing ⚠️
...solve/nonlinear/descent_strategies/NonlinearCG.cpp 71.69% 15 Missing ⚠️
src/polysolve/linear/HypreSolver.cpp 20.00% 8 Missing ⚠️
src/polysolve/nonlinear/Solver.cpp 58.82% 7 Missing ⚠️
.../polysolve/nonlinear/descent_strategies/Newton.hpp 0.00% 3 Missing ⚠️
src/polysolve/linear/CPUHybridSolver.hpp 50.00% 1 Missing ⚠️
src/polysolve/linear/Solver.hpp 50.00% 1 Missing ⚠️
src/polysolve/nonlinear/Problem.hpp 75.00% 1 Missing ⚠️
...rc/polysolve/nonlinear/descent_strategies/BFGS.cpp 80.00% 1 Missing ⚠️

❌ Your patch check has failed because the patch coverage (51.13%) is below the target coverage (75.00%). You can increase the patch coverage or adjust the target coverage.

Additional details and impacted files
@@            Coverage Diff             @@
##             main     #120      +/-   ##
==========================================
- Coverage   79.45%   75.95%   -3.50%     
==========================================
  Files          55       57       +2     
  Lines        2940     3186     +246     
  Branches      377      409      +32     
==========================================
+ Hits         2336     2420      +84     
- Misses        604      766     +162     
Flag Coverage Δ
polysolve 75.95% <51.13%> (-3.50%) ⬇️

Flags with carried forward coverage won't be shown. Click here to find out more.

☔ View full report in Codecov by Harness.
📢 Have feedback on the report? Share it here.

🚀 New features to boost your workflow:
  • ❄️ Test Analytics: Detect flaky tests, report on failures, and find test suite problems.

Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment

Labels

None yet

Projects

None yet

Development

Successfully merging this pull request may close these issues.

2 participants