`reuse_preconditioner` destroys the matrix-free operator on a retained SNES: `-snes_mf_operator` silently becomes a Newton solve on the preconditioning matrix, `-snes_mf` can crash

[Lewin-Shen]

With PetscNonlinearSolver::reuse_preconditioner() on, the SNES is kept alive between solves. On every solve,
PetscNonlinearSolver<T>::solve() registers the matrices it is given as the SNES operators,
SNESSetJacobian (_snes, jac->mat(), pre->mat(), libmesh_petsc_snes_jacobian, this) (src/solvers/petsc_nonlinear_solver.C:965-968
at devel cf92df06d9). The single-matrix solve(), which NonlinearImplicitSystem::solve() calls unless an operator matrix was
set with set_operator_matrix(), passes the assembled preconditioning matrix in both slots (:900-909), so that matrix is
registered as both operators, as it was before #4537 (:951-954 at 90766057ef). On a retained SNES the SNESSetUp() that
follows returns at its early exit (snes.c: if (snes->setupcalled) PetscFunctionReturn(PETSC_SUCCESS);), so the matrix-free
operator that -snes_mf_operator or -snes_mf asked for, and that the first solve's SNESSetUp() created, is never re-created.
The registration drops the SNES's reference to it and the options step that follows drops the KSP's, so it is destroyed.

Consequences, observed at 90766057ef with PETSc 3.25.4 at 1 and 2 ranks, and reproduced by the pull request's unit test at
devel cf92df06d9:

• -snes_mf_operator (MOOSE's PJFNK): from the second solve on, the Krylov method is applied to the assembled preconditioning
  matrix. The run completes and gives a different answer with no message about the operator: the solves that fail read as a
  time-step problem, and -snes_view still prints Jacobian is applied matrix-free with differencing, because that line reports
  the flag and not the installed matrix. In one MOOSE simulation the time step collapsed from 20 to 0.625, the run stopped at
  simulated time 25 where the same run without reuse reaches 120, and a reported integral differed by 9 %; the exit status was 0.
• -snes_mf: the Jacobian callback PETSc installed (MatMFFDComputeJacobian) is lost as well, because init() attaches a new
  DMLIBMESH on every solve, DMSetUp_libMesh registers its own DMSNES callbacks (src/solvers/petscdmlibmeshimpl.C:905-906),
  and SNESSetDM copies the old DMSNES only onto a DM that has none (snes.c: if (snes->dm->dmsnes && !dm->dmsnes)). With only
  a residual registered, the DM's own Jacobian callback then receives PETSc's matrix-free matrix, PetscMatrixBase::get_context()
  returns null for it, and the process dies with a segmentation fault in SparseMatrix::attach_dof_map. With a Jacobian
  registered as well, the second solve becomes a Newton solve on the assembled matrix instead.

NEWTON solves (assembled operator) are not affected. No test in libMesh or MOOSE runs two solves with reuse on and a
matrix-free solve type, and the first solve is correct by construction, which is why every existing test passes.

Reproducer: two solves of one NonlinearImplicitSystem with reuse preconditioner set through the EquationSystems
parameters, under -snes_mf_operator (or -snes_mf), inspecting SNESGetJacobian after each solve: solve 1 has a mffd
operator, solve 2 has the assembled matrix (or crashes). The pull request below adds this as a unit test.

Affected versions: demonstrated at libMesh 90766057ef and at devel cf92df06d9 (2026-09-24) with PETSc 3.25.4. At
devel the code path is unchanged: #4537, merged on 2026-09-24, moved the registration into the new two-matrix solve() and left the
single-matrix solve() registering the preconditioning matrix as both operators. PETSc main keeps the early exit and the
matrix-free placement. Suspected, pending a history read: every libMesh since the reuse feature landed in 2022 (#3208, #3209),
because the registration this depends on predates the feature.

Proposed fix: a pull request follows: keep the retained matrix-free operator (about 50 lines of code in solve()).