From 9b9ba71cb8eb705e9ed0192d2c1ea78e1ebfd79a Mon Sep 17 00:00:00 2001 From: Nuno Nobre Date: Wed, 25 Feb 2026 12:06:36 +0000 Subject: [PATCH 1/4] 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 4cc0ccd6d4..2f56b55de4 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 ec18fa24382eb396ddc91bf16739fa7fec9c501b Mon Sep 17 00:00:00 2001 From: Nuno Nobre Date: Fri, 7 Feb 2025 19:24:03 +0000 Subject: [PATCH 2/4] 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 9a7ded89a6..39d19fa0a8 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 2f56b55de4..74d038c835 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 1261cf01cfd91763e9713eb1c85a47eeb172fd43 Mon Sep 17 00:00:00 2001 From: Nuno Nobre Date: Wed, 25 Feb 2026 15:32:21 +0000 Subject: [PATCH 3/4] Prevent nested FSP as I dunno how to get matrix blocks --- src/numerics/petsc_preconditioner.C | 15 ++++++++++++++- 1 file changed, 14 insertions(+), 1 deletion(-) diff --git a/src/numerics/petsc_preconditioner.C b/src/numerics/petsc_preconditioner.C index 39d19fa0a8..85a20337cd 100644 --- a/src/numerics/petsc_preconditioner.C +++ b/src/numerics/petsc_preconditioner.C @@ -173,6 +173,12 @@ void PetscPreconditioner::set_petsc_aux_data(PC & pc, System & sys, const uns PCType pc_type = nullptr; LibmeshPetscCallA(comm, PCGetType(pc, &pc_type)); +#if !PETSC_VERSION_LESS_THAN(3, 20, 0) + // Get the nesting level of this PC's KSP + PetscInt level; + LibmeshPetscCallA(comm, PCGetKSPNestLevel(pc, &level)); +#endif + // Check if hypre ams/ads, otherwise we quit with nothing to do if (pc_type && std::string(pc_type) == PCHYPRE) { @@ -205,7 +211,8 @@ 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) +#if !PETSC_VERSION_LESS_THAN(3, 20, 0) + 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,17 +226,23 @@ 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)); + LibmeshPetscCallA(comm, PCSetKSPNestLevel(subpc, level + 1)); set_petsc_aux_data(subpc, sys, s); } // Free the array of sub KSP contexts LibmeshPetscCallA(comm, PetscFree(subksps)); } +#endif } From 71c9d1985124a54fd0d1bcbe36aff75cb3a75f0e Mon Sep 17 00:00:00 2001 From: Nuno Nobre Date: Sat, 11 Apr 2026 13:49:16 +0100 Subject: [PATCH 4/4] Disable aux data with PETSc DM [do not merge] --- src/solvers/petsc_linear_solver.C | 7 +++++-- 1 file changed, 5 insertions(+), 2 deletions(-) diff --git a/src/solvers/petsc_linear_solver.C b/src/solvers/petsc_linear_solver.C index ae8c813acd..ad7a8c57c3 100644 --- a/src/solvers/petsc_linear_solver.C +++ b/src/solvers/petsc_linear_solver.C @@ -223,8 +223,11 @@ template void PetscLinearSolver::init_systems (const System & sys) { - petsc_auto_fieldsplit(this->pc(), sys); - PetscPreconditioner::set_petsc_aux_data(_pc, const_cast(sys)); + if(!libMesh::on_command_line("--use_petsc_dm")) + { + petsc_auto_fieldsplit(this->pc(), sys); + PetscPreconditioner::set_petsc_aux_data(_pc, const_cast(sys)); + } }