diff --git a/src/numerics/petsc_preconditioner.C b/src/numerics/petsc_preconditioner.C index 2c700d3821..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) { @@ -183,10 +187,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 +199,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 +209,38 @@ 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 && !level) + { + // 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)); + + // 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); + } + + // 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 e2a09985ec..89cae9fc80 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) || @@ -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) {