From fa0c1903642c0dea1ea358d5d725a35d175e8343 Mon Sep 17 00:00:00 2001 From: Nuno Nobre Date: Wed, 25 Feb 2026 12:06:36 +0000 Subject: [PATCH 1/3] Rename PC->Pc to avoid hiding typename (now mimics jac -> Jac) --- src/solvers/petsc_nonlinear_solver.C | 22 +++++++++++----------- 1 file changed, 11 insertions(+), 11 deletions(-) diff --git a/src/solvers/petsc_nonlinear_solver.C b/src/solvers/petsc_nonlinear_solver.C index e2a09985ec..04aa98e328 100644 --- a/src/solvers/petsc_nonlinear_solver.C +++ b/src/solvers/petsc_nonlinear_solver.C @@ -432,7 +432,7 @@ extern "C" NonlinearImplicitSystem & sys = solver->system(); - PetscMatrixBase * const PC = pc ? PetscMatrixBase::get_context(pc, sys.comm()) : nullptr; + PetscMatrixBase * const Pc = pc ? PetscMatrixBase::get_context(pc, sys.comm()) : nullptr; PetscMatrixBase * Jac = jac ? PetscMatrixBase::get_context(jac, sys.comm()) : nullptr; PetscVector & X_sys = *cast_ptr *>(sys.solution.get()); PetscVector X_global(x, sys.comm()); @@ -462,13 +462,13 @@ extern "C" Jac->close(); if (pc && (p_is_shell == PETSC_TRUE)) - PC->close(); + Pc->close(); PetscFunctionReturn(LIBMESH_PETSC_SUCCESS); } // Set the dof maps - PC->attach_dof_map(sys.get_dof_map()); + Pc->attach_dof_map(sys.get_dof_map()); Jac->attach_dof_map(sys.get_dof_map()); // Use the systems update() to get a good local version of the parallel solution @@ -484,29 +484,29 @@ extern "C" sys.get_dof_map().enforce_constraints_exactly(sys, sys.current_local_solution.get()); if (solver->_zero_out_jacobian) - PC->zero(); + Pc->zero(); if (solver->jacobian != nullptr) - solver->jacobian(*sys.current_local_solution.get(), *PC, sys); + solver->jacobian(*sys.current_local_solution.get(), *Pc, sys); else if (solver->jacobian_object != nullptr) - solver->jacobian_object->jacobian(*sys.current_local_solution.get(), *PC, sys); + solver->jacobian_object->jacobian(*sys.current_local_solution.get(), *Pc, sys); else if (solver->matvec != nullptr) - solver->matvec(*sys.current_local_solution.get(), nullptr, PC, sys); + solver->matvec(*sys.current_local_solution.get(), nullptr, Pc, sys); else libmesh_error_msg("Error! Unable to compute residual and/or Jacobian!"); - PC->close(); + Pc->close(); if (solver->_exact_constraint_enforcement) { - sys.get_dof_map().enforce_constraints_on_jacobian(sys, PC); - PC->close(); + sys.get_dof_map().enforce_constraints_on_jacobian(sys, Pc); + Pc->close(); } - if (Jac != PC) + if (Jac != Pc) { // Assume that shells know what they're doing libmesh_assert(!solver->_exact_constraint_enforcement || (j_is_mffd == PETSC_TRUE) || From 89899fe7f417cf7c8e2eba5b128a00c0f01e41ef Mon Sep 17 00:00:00 2001 From: Nuno Nobre Date: Fri, 7 Feb 2025 19:24:03 +0000 Subject: [PATCH 2/3] Allow integrating hypre AMS/ADS in FSPs --- src/numerics/petsc_preconditioner.C | 33 +++++++++++++++++++++------- src/solvers/petsc_nonlinear_solver.C | 10 +++++---- 2 files changed, 31 insertions(+), 12 deletions(-) diff --git a/src/numerics/petsc_preconditioner.C b/src/numerics/petsc_preconditioner.C index 2c700d3821..9ccef4840e 100644 --- a/src/numerics/petsc_preconditioner.C +++ b/src/numerics/petsc_preconditioner.C @@ -183,10 +183,6 @@ void PetscPreconditioner::set_petsc_aux_data(PC & pc, System & sys, const uns // If not ams/ads, we quit with nothing to do if (std::string(hypre_type) == "ams") { - // If multiple variables, we error out as senseless - libmesh_error_msg_if(sys.n_vars() > 1, - "Error applying hypre AMS to a system with multiple " - "variables"); // If not a 1st order Nédélec or a 2d 1st order Raviart-Thomas system, we // error out as we do not support anything else at the moment libmesh_error_msg_if(sys.variable(v).type() != FEType(1, NEDELEC_ONE) && @@ -199,10 +195,6 @@ void PetscPreconditioner::set_petsc_aux_data(PC & pc, System & sys, const uns } else if (std::string(hypre_type) == "ads") { - // If multiple variables, we error out as senseless - libmesh_error_msg_if(sys.n_vars() > 1, - "Error applying hypre ADS to a system with multiple " - "variables"); // If not a 3d 1st order Raviart-Thomas system, we error out as we do not // support anything else at the moment libmesh_error_msg_if(sys.variable(v).type() != FEType(1, RAVIART_THOMAS) || @@ -213,6 +205,31 @@ void PetscPreconditioner::set_petsc_aux_data(PC & pc, System & sys, const uns set_hypre_ads_data(pc, sys, v); } } + else if (pc_type && std::string(pc_type) == PCFIELDSPLIT) + { + // Annoyingly, if using Schur complement preconditioning, we need to call PCSetUp() + auto pmatrix = cast_ptr *>(sys.request_matrix("System Matrix")); + pmatrix->close(); + Mat mat = pmatrix->mat(); + LibmeshPetscCallA(comm, PCSetOperators(pc, mat, mat)); + LibmeshPetscCallA(comm, PCSetUp(pc)); + + // Get sub KSP contexts for all splits + PetscInt n_splits; + KSP * subksps; + LibmeshPetscCallA(comm, PCFieldSplitGetSubKSP(pc, &n_splits, &subksps)); + + // Get the sub PC context for each split and recursively call this function + for (auto s : make_range(n_splits)) + { + PC subpc; + LibmeshPetscCallA(comm, KSPGetPC(subksps[s], &subpc)); + set_petsc_aux_data(subpc, sys, s); + } + + // Free the array of sub KSP contexts + LibmeshPetscCallA(comm, PetscFree(subksps)); + } } diff --git a/src/solvers/petsc_nonlinear_solver.C b/src/solvers/petsc_nonlinear_solver.C index 04aa98e328..89cae9fc80 100644 --- a/src/solvers/petsc_nonlinear_solver.C +++ b/src/solvers/petsc_nonlinear_solver.C @@ -514,6 +514,12 @@ extern "C" Jac->close(); } + KSP kspc; + PC pcc; + LibmeshPetscCall2(sys.comm(), SNESGetKSP(snes, &kspc)); + LibmeshPetscCall2(sys.comm(), KSPGetPC(kspc, &pcc)); + PetscPreconditioner::set_petsc_aux_data(pcc, sys); + PetscFunctionReturn(LIBMESH_PETSC_SUCCESS); } @@ -985,10 +991,6 @@ PetscNonlinearSolver::solve (SparseMatrix & pre_in, // System Preconditi #endif LibmeshPetscCall(SNESSetFromOptions(_snes)); - PC pc; - LibmeshPetscCall(KSPGetPC(ksp, &pc)); - PetscPreconditioner::set_petsc_aux_data(pc, this->system()); - #if defined(LIBMESH_HAVE_PETSC_HYPRE) && PETSC_VERSION_LESS_THAN(3, 23, 0) && \ !PETSC_VERSION_LESS_THAN(3, 12, 0) && defined(PETSC_HAVE_HYPRE_DEVICE) { From d01f48cdc7b6bbf1fde7bd55a29f633c761af628 Mon Sep 17 00:00:00 2001 From: Nuno Nobre Date: Wed, 25 Feb 2026 15:32:21 +0000 Subject: [PATCH 3/3] Prevent nested FSP as I dunno how to get matrix blocks --- src/numerics/petsc_preconditioner.C | 13 ++++++++++++- 1 file changed, 12 insertions(+), 1 deletion(-) diff --git a/src/numerics/petsc_preconditioner.C b/src/numerics/petsc_preconditioner.C index 9ccef4840e..b0c3a8e078 100644 --- a/src/numerics/petsc_preconditioner.C +++ b/src/numerics/petsc_preconditioner.C @@ -173,6 +173,10 @@ void PetscPreconditioner::set_petsc_aux_data(PC & pc, System & sys, const uns PCType pc_type = nullptr; LibmeshPetscCallA(comm, PCGetType(pc, &pc_type)); + // Get the nesting level of this PC's KSP + PetscInt level; + LibmeshPetscCallA(comm, PCGetKSPNestLevel(pc, &level)); + // Check if hypre ams/ads, otherwise we quit with nothing to do if (pc_type && std::string(pc_type) == PCHYPRE) { @@ -205,7 +209,7 @@ void PetscPreconditioner::set_petsc_aux_data(PC & pc, System & sys, const uns set_hypre_ads_data(pc, sys, v); } } - else if (pc_type && std::string(pc_type) == PCFIELDSPLIT) + else if (pc_type && std::string(pc_type) == PCFIELDSPLIT && !level) { // Annoyingly, if using Schur complement preconditioning, we need to call PCSetUp() auto pmatrix = cast_ptr *>(sys.request_matrix("System Matrix")); @@ -219,11 +223,18 @@ void PetscPreconditioner::set_petsc_aux_data(PC & pc, System & sys, const uns KSP * subksps; LibmeshPetscCallA(comm, PCFieldSplitGetSubKSP(pc, &n_splits, &subksps)); + // We assume a one-to-one map between splits and variables + if (sys.n_vars() != n_splits) + return; + // Get the sub PC context for each split and recursively call this function for (auto s : make_range(n_splits)) { PC subpc; LibmeshPetscCallA(comm, KSPGetPC(subksps[s], &subpc)); +#if PETSC_VERSION_LESS_THAN(3, 25, 0) + LibmeshPetscCallA(comm, PCSetKSPNestLevel(subpc, level + 1)); +#endif set_petsc_aux_data(subpc, sys, s); }