PetscSupport.C

Go to the documentation of this file.

//* This file is part of the MOOSE framework
//* https://mooseframework.inl.gov
//*
//* All rights reserved, see COPYRIGHT for full restrictions
//* https://github.com/idaholab/moose/blob/master/COPYRIGHT
//*
//* Licensed under LGPL 2.1, please see LICENSE for details
//* https://www.gnu.org/licenses/lgpl-2.1.html
 
#include "PetscSupport.h"
 
// MOOSE includes
#include "MooseApp.h"
#include "FEProblem.h"
#include "DisplacedProblem.h"
#include "NonlinearSystem.h"
#include "LinearSystem.h"
#include "AuxiliarySystem.h"
#include "DisplacedProblem.h"
#include "PenetrationLocator.h"
#include "NearestNodeLocator.h"
#include "MooseTypes.h"
#include "MooseUtils.h"
#include "CommandLine.h"
#include "Console.h"
#include "MultiMooseEnum.h"
#include "Conversion.h"
#include "Executioner.h"
#include "MooseMesh.h"
#include "ComputeLineSearchObjectWrapper.h"
#include "Convergence.h"
#include "ParallelParamObject.h"
 
#include "libmesh/equation_systems.h"
#include "libmesh/linear_implicit_system.h"
#include "libmesh/nonlinear_implicit_system.h"
#include "libmesh/petsc_linear_solver.h"
#include "libmesh/petsc_matrix.h"
#include "libmesh/petsc_nonlinear_solver.h"
#include "libmesh/petsc_preconditioner.h"
#include "libmesh/petsc_vector.h"
#include "libmesh/sparse_matrix.h"
#include "libmesh/petsc_solver_exception.h"
#include "libmesh/simple_range.h"
 
// PETSc includes
#include <petsc.h>
#include <petscsnes.h>
#include <petscksp.h>
#include <petscmat.h>
#include <petscis.h>
#include <petscdm.h>
#include <petscoptions.h>
 
// PetscDMMoose include
#include "PetscDMMoose.h"
 
// Standard includes
#include <ostream>
#include <cctype>
#include <fstream>
#include <optional>
#include <string>
 
using namespace libMesh;
 
void
MooseVecView(NumericVector<Number> & vector)
{
  PetscVector<Number> & petsc_vec = static_cast<PetscVector<Number> &>(vector);
  LibmeshPetscCallA(vector.comm().get(), VecView(petsc_vec.vec(), 0));
}
 
void
MooseMatView(SparseMatrix<Number> & mat)
{
  PetscMatrixBase<Number> & petsc_mat = static_cast<PetscMatrix<Number> &>(mat);
  LibmeshPetscCallA(mat.comm().get(), MatView(petsc_mat.mat(), 0));
}
 
void
MooseVecView(const NumericVector<Number> & vector)
{
  PetscVector<Number> & petsc_vec =
      static_cast<PetscVector<Number> &>(const_cast<NumericVector<Number> &>(vector));
  LibmeshPetscCallA(vector.comm().get(), VecView(petsc_vec.vec(), 0));
}
 
void
MooseMatView(const SparseMatrix<Number> & mat)
{
  PetscMatrixBase<Number> & petsc_mat =
      static_cast<PetscMatrix<Number> &>(const_cast<SparseMatrix<Number> &>(mat));
  LibmeshPetscCallA(mat.comm().get(), MatView(petsc_mat.mat(), 0));
}
 
namespace Moose
{
namespace PetscSupport
{
 
PetscOptionsScope::PetscOptionsScope(FEProblemBase & problem) : _problem(problem), _pushed(false)
{
#if !PETSC_RELEASE_LESS_THAN(3, 12, 0)
  if (!_problem.getMooseApp().isUltimateMaster())
  {
    LibmeshPetscCallA(_problem.comm().get(), PetscOptionsPush(_problem.petscOptionsDatabase()));
    _pushed = true;
  }
#endif
}
 
PetscOptionsScope::~PetscOptionsScope()
{
#if !PETSC_RELEASE_LESS_THAN(3, 12, 0)
  if (_pushed)
    PetscCallAbort(_problem.comm().get(), PetscOptionsPop());
#endif
}
 
namespace
{
 
void
applySystemVectorTypeOptions(FEProblemBase & problem, libMesh::System & lm_sys)
{
  for (auto & [_, vec] : as_range(lm_sys.vectors_begin(), lm_sys.vectors_end()))
  {
    auto * const petsc_vec = cast_ptr<PetscVector<Number> *>(vec.get());
    LibmeshPetscCallA(problem.comm().get(), VecSetFromOptions(petsc_vec->vec()));
  }
 
  // The solution vectors aren't included in the system vectors storage.
  auto * petsc_vec = cast_ptr<PetscVector<Number> *>(lm_sys.solution.get());
  LibmeshPetscCallA(problem.comm().get(), VecSetFromOptions(petsc_vec->vec()));
  petsc_vec = cast_ptr<PetscVector<Number> *>(lm_sys.current_local_solution.get());
  LibmeshPetscCallA(problem.comm().get(), VecSetFromOptions(petsc_vec->vec()));
}
 
void
applyVectorTypeOptions(FEProblemBase & problem)
{
  for (const auto sys_index : make_range(problem.numSolverSystems()))
    applySystemVectorTypeOptions(problem, problem.getSolverSystem(sys_index).system());
 
  applySystemVectorTypeOptions(problem, problem.getAuxiliarySystem().system());
}
 
bool
petscOptionsHasName(::PetscOptions options,
                    const std::string & name,
                    const std::string & prefix = "")
{
  PetscBool found = PETSC_FALSE;
  const char * const prefix_ptr = prefix.empty() ? nullptr : prefix.c_str();
  LibmeshPetscCallA(PETSC_COMM_WORLD,
                    PetscOptionsHasName(options, prefix_ptr, name.c_str(), &found));
  return found;
}
 
bool
hasMatrixFreeSolveType(const FEProblemBase & problem)
{
  for (const auto sys_index : make_range(problem.numSolverSystems()))
    if (const auto solve_type = problem.solverParams(sys_index)._type;
        solve_type == Moose::ST_JFNK || solve_type == Moose::ST_PJFNK)
      return true;
 
  return false;
}
 
bool
mightBeMatTypeOption(const std::string & name)
{
  static constexpr std::string_view mat_type_suffix = "mat_type";
  return name.size() >= mat_type_suffix.size() &&
         std::equal(mat_type_suffix.rbegin(),
                    mat_type_suffix.rend(),
                    name.rbegin(),
                    [](const char left, const char right)
                    // tolower requires representability by unsigned char
                    {
                      return static_cast<int>(left) ==
                             std::tolower(libMesh::cast_int<unsigned char>(right));
                    });
}
 
void
errorOnUnprefixedMatTypeOption(::PetscOptions options, FEProblemBase & problem)
{
  if (!petscOptionsHasName(options, "-mat_type"))
    return;
 
  std::string error_string =
      "Setting option '-mat_type' is not supported without a solver-system prefix. Use an option "
      "such as '-" +
      problem.getSolverSystem(0).name() + "_mat_type' for assembled libMesh matrices.";
  if (hasMatrixFreeSolveType(problem))
    error_string +=
        " Attempting to change the matrix "
        "type for the MFFD matrix type used to represent the Jacobian for (P)JFNK solve "
        "types is not supported.";
  mooseError(error_string);
}
 
// Allow iterating over all systems or allowing caller to specify a specific system for which to
// apply matrix type options
void
applyMatrixTypeOptions(FEProblemBase & problem,
                       const std::optional<std::size_t> system_index = std::nullopt)
{
  const auto begin = system_index.value_or(0);
  const auto end = system_index ? begin + 1 : problem.numSolverSystems();
 
  for (const auto sys_index : make_range(begin, end))
  {
    auto & solver_system = problem.getSolverSystem(sys_index);
    auto & lm_sys = solver_system.system();
 
    // Even in matrix-free modes the libMesh matrix wrappers can exist before the PETSc Mat does.
    if (problem.solverParams(sys_index)._type == Moose::ST_JFNK)
      continue;
 
    for (auto & [_, mat] : as_range(lm_sys.matrices_begin(), lm_sys.matrices_end()))
      if (auto * const petsc_mat = dynamic_cast<PetscMatrixBase<Number> *>(mat.get()); petsc_mat)
      {
        LibmeshPetscCallA(
            problem.comm().get(),
            MatSetOptionsPrefix(petsc_mat->mat(), (solver_system.name() + "_").c_str()));
        LibmeshPetscCallA(problem.comm().get(), MatSetFromOptions(petsc_mat->mat()));
      }
  }
}
 
} // namespace
 
std::string
stringify(const LineSearchType & t)
{
  switch (t)
  {
    case LS_BASIC:
      return "basic";
    case LS_DEFAULT:
      return "default";
    case LS_NONE:
      return "none";
    case LS_SHELL:
      return "shell";
    case LS_L2:
      return "l2";
    case LS_BT:
      return "bt";
    case LS_CP:
      return "cp";
    case LS_CONTACT:
      return "contact";
    case LS_PROJECT:
      return "project";
    case LS_INVALID:
      mooseError("Invalid LineSearchType");
  }
  return "";
}
 
std::string
stringify(const MffdType & t)
{
  switch (t)
  {
    case MFFD_WP:
      return "wp";
    case MFFD_DS:
      return "ds";
    case MFFD_INVALID:
      mooseError("Invalid MffdType");
  }
  return "";
}
 
void
setSolverOptions(const SolverParams & solver_params, const MultiMooseEnum & dont_add_these_options)
{
  const auto prefix_with_dash = '-' + solver_params._prefix;
  // set PETSc options implied by a solve type
  switch (solver_params._type)
  {
    case Moose::ST_PJFNK:
      setSinglePetscOptionIfAppropriate(dont_add_these_options,
                                        prefix_with_dash + "snes_mf_operator");
      setSinglePetscOptionIfAppropriate(dont_add_these_options,
                                        prefix_with_dash + "mat_mffd_type",
                                        stringify(solver_params._mffd_type));
      break;
 
    case Moose::ST_JFNK:
      setSinglePetscOptionIfAppropriate(dont_add_these_options, prefix_with_dash + "snes_mf");
      setSinglePetscOptionIfAppropriate(dont_add_these_options,
                                        prefix_with_dash + "mat_mffd_type",
                                        stringify(solver_params._mffd_type));
      break;
 
    case Moose::ST_NEWTON:
      break;
 
    case Moose::ST_FD:
      setSinglePetscOptionIfAppropriate(dont_add_these_options, prefix_with_dash + "snes_fd");
      break;
 
    case Moose::ST_LINEAR:
      setSinglePetscOptionIfAppropriate(
          dont_add_these_options, prefix_with_dash + "snes_type", "ksponly");
      setSinglePetscOptionIfAppropriate(dont_add_these_options,
                                        prefix_with_dash + "snes_monitor_cancel");
      break;
  }
 
  Moose::LineSearchType ls_type = solver_params._line_search;
  if (ls_type == Moose::LS_NONE)
    ls_type = Moose::LS_BASIC;
 
  if (ls_type != Moose::LS_DEFAULT && ls_type != Moose::LS_CONTACT && ls_type != Moose::LS_PROJECT)
    setSinglePetscOptionIfAppropriate(
        dont_add_these_options, prefix_with_dash + "snes_linesearch_type", stringify(ls_type));
}
 
void
addPetscOptionsFromCommandline(FEProblemBase * const problem)
{
  // commandline options always win
  // the options from a user commandline will overwrite the existing ones if any conflicts
  int argc;
  char ** args;
 
  LibmeshPetscCallA(PETSC_COMM_WORLD, PetscGetArgs(&argc, &args));
  std::vector<const char *> cl_args(args + 1, args + argc);
  const auto cl_argc = libMesh::cast_int<int>(cl_args.size());
 
  ::PetscOptions command_line_options;
  LibmeshPetscCallA(PETSC_COMM_WORLD, PetscOptionsCreate(&command_line_options));
  LibmeshPetscCallA(PETSC_COMM_WORLD,
                    PetscOptionsInsertArgs(command_line_options, cl_argc, cl_args.data()));
  LibmeshPetscCallA(PETSC_COMM_WORLD,
                    PetscOptionsInsertArgs(LIBMESH_PETSC_NULLPTR, cl_argc, cl_args.data()));
 
  if (!problem)
  {
    LibmeshPetscCallA(PETSC_COMM_WORLD, PetscOptionsDestroy(&command_line_options));
    return;
  }
 
  errorOnUnprefixedMatTypeOption(command_line_options, *problem);
 
  // Some vector/matrix-type options may have been consumed before the PETSc database rebuild.
  // Replay only the command-line-controlled applications so input-file options handled through
  // setSinglePetscOption() do not pay the cost twice.
  const bool have_vec_type = petscOptionsHasName(command_line_options, "-vec_type");
  bool have_mat_type = false;
 
  for (const auto sys_index : make_range(problem->numSolverSystems()))
  {
    have_mat_type = petscOptionsHasName(
        command_line_options, "-mat_type", problem->getSolverSystem(sys_index).name() + "_");
    if (have_mat_type)
      break;
  }
 
  if (have_vec_type)
    applyVectorTypeOptions(*problem);
  if (have_mat_type)
    applyMatrixTypeOptions(*problem);
 
  LibmeshPetscCallA(PETSC_COMM_WORLD, PetscOptionsDestroy(&command_line_options));
}
 
void
petscSetOptionsHelper(const PetscOptions & po, FEProblemBase * const problem)
{
  // Add any additional options specified in the input file
  for (const auto & flag : po.flags)
    // Need to use name method here to pass a str instead of an EnumItem because
    // we don't care if the id attributes match
    if (!po.dont_add_these_options.contains(flag.name()) ||
        po.user_set_options.contains(flag.name()))
      setSinglePetscOption(flag.rawName().c_str());
 
  // Add option pairs
  for (auto & option : po.pairs)
    if (!po.dont_add_these_options.contains(option.first) ||
        po.user_set_options.contains(option.first))
      setSinglePetscOption(option.first, option.second, problem);
 
  addPetscOptionsFromCommandline(problem);
}
 
void
petscSetOptions(const PetscOptions & po,
                const SolverParams & solver_params,
                FEProblemBase * const problem)
{
  PetscCallAbort(PETSC_COMM_WORLD, PetscOptionsClear(LIBMESH_PETSC_NULLPTR));
  setSolverOptions(solver_params, po.dont_add_these_options);
  petscSetOptionsHelper(po, problem);
}
 
void
petscSetOptions(const PetscOptions & po,
                const std::vector<SolverParams> & solver_params_vec,
                FEProblemBase * const problem)
{
  PetscCallAbort(PETSC_COMM_WORLD, PetscOptionsClear(LIBMESH_PETSC_NULLPTR));
  for (const auto & solver_params : solver_params_vec)
    setSolverOptions(solver_params, po.dont_add_these_options);
  petscSetOptionsHelper(po, problem);
}
 
PetscErrorCode
petscSetupOutput(CommandLine * cmd_line)
{
  PetscFunctionBegin;
  char code[10] = {45, 45, 109, 111, 111, 115, 101};
  const std::vector<std::string> argv = cmd_line->getArguments();
  for (const auto & arg : argv)
  {
    if (arg.compare(code) == 0)
    {
      Console::petscSetupOutput();
      break;
    }
  }
  PetscFunctionReturn(PETSC_SUCCESS);
}
 
PetscErrorCode
petscNonlinearConverged(SNES /*snes*/,
                        PetscInt it,
                        PetscReal /*xnorm*/,
                        PetscReal /*snorm*/,
                        PetscReal /*fnorm*/,
                        SNESConvergedReason * reason,
                        void * ctx)
{
  PetscFunctionBegin;
  FEProblemBase & problem = *static_cast<FEProblemBase *>(ctx);
 
  // execute objects that may be used in convergence check
  problem.execute(EXEC_NONLINEAR_CONVERGENCE);
 
  // perform the convergence check
  Convergence::MooseConvergenceStatus status;
  if (problem.getFailNextNonlinearConvergenceCheck())
  {
    status = Convergence::MooseConvergenceStatus::DIVERGED;
    problem.resetFailNextNonlinearConvergenceCheck();
  }
  else
  {
    auto & convergence = problem.getConvergence(
        problem.getNonlinearConvergenceNames()[problem.currentNonlinearSystem().number()]);
    status = convergence.checkConvergence(it);
  }
 
  // convert convergence status to PETSc converged reason
  switch (status)
  {
    case Convergence::MooseConvergenceStatus::ITERATING:
      *reason = SNES_CONVERGED_ITERATING;
      break;
 
    case Convergence::MooseConvergenceStatus::CONVERGED:
      *reason = SNES_CONVERGED_FNORM_ABS;
      break;
 
    case Convergence::MooseConvergenceStatus::DIVERGED:
      *reason = SNES_DIVERGED_DTOL;
      break;
  }
 
  PetscFunctionReturn(PETSC_SUCCESS);
}
 
PetscErrorCode
petscLinearConverged(
    KSP /*ksp*/, PetscInt it, PetscReal /*norm*/, KSPConvergedReason * reason, void * ctx)
{
  PetscFunctionBegin;
  FEProblemBase & problem = *static_cast<FEProblemBase *>(ctx);
 
  // execute objects that may be used in convergence check
  // Right now, setting objects to execute on this flag would be ignored except in the
  // linear-system-only use case.
  problem.execute(EXEC_LINEAR_CONVERGENCE);
 
  // perform the convergence check
  Convergence::MooseConvergenceStatus status;
  if (problem.getFailNextSystemConvergenceCheck())
  {
    status = Convergence::MooseConvergenceStatus::DIVERGED;
    problem.resetFailNextSystemConvergenceCheck();
  }
  else
  {
    auto & convergence = problem.getConvergence(
        problem.getLinearConvergenceNames()[problem.currentLinearSystem().number()]);
    status = convergence.checkConvergence(it);
  }
 
  // convert convergence status to PETSc converged reason
  switch (status)
  {
    case Convergence::MooseConvergenceStatus::ITERATING:
      *reason = KSP_CONVERGED_ITERATING;
      break;
 
      // TODO: find a KSP code that works better for this case
    case Convergence::MooseConvergenceStatus::CONVERGED:
#if PETSC_VERSION_LESS_THAN(3, 24, 0)
      *reason = KSP_CONVERGED_RTOL_NORMAL;
#else
      *reason = KSP_CONVERGED_RTOL_NORMAL_EQUATIONS;
#endif
      break;
 
    case Convergence::MooseConvergenceStatus::DIVERGED:
      *reason = KSP_DIVERGED_DTOL;
      break;
  }
 
  PetscFunctionReturn(PETSC_SUCCESS);
}
 
PCSide
getPetscPCSide(Moose::PCSideType pcs)
{
  switch (pcs)
  {
    case Moose::PCS_LEFT:
      return PC_LEFT;
    case Moose::PCS_RIGHT:
      return PC_RIGHT;
    case Moose::PCS_SYMMETRIC:
      return PC_SYMMETRIC;
    default:
      mooseError("Unknown PC side requested.");
      break;
  }
}
 
KSPNormType
getPetscKSPNormType(Moose::MooseKSPNormType kspnorm)
{
  switch (kspnorm)
  {
    case Moose::KSPN_NONE:
      return KSP_NORM_NONE;
    case Moose::KSPN_PRECONDITIONED:
      return KSP_NORM_PRECONDITIONED;
    case Moose::KSPN_UNPRECONDITIONED:
      return KSP_NORM_UNPRECONDITIONED;
    case Moose::KSPN_NATURAL:
      return KSP_NORM_NATURAL;
    case Moose::KSPN_DEFAULT:
      return KSP_NORM_DEFAULT;
    default:
      mooseError("Unknown KSP norm type requested.");
      break;
  }
}
 
void
petscSetDefaultKSPNormType(FEProblemBase & problem, KSP ksp)
{
  for (const auto i : make_range(problem.numSolverSystems()))
  {
    SolverSystem & sys = problem.getSolverSystem(i);
    LibmeshPetscCallA(problem.comm().get(),
                      KSPSetNormType(ksp, getPetscKSPNormType(sys.getMooseKSPNormType())));
  }
}
 
void
petscSetDefaultPCSide(FEProblemBase & problem, KSP ksp)
{
  for (const auto i : make_range(problem.numSolverSystems()))
  {
    SolverSystem & sys = problem.getSolverSystem(i);
 
    // PETSc 3.2.x+
    if (sys.getPCSide() != Moose::PCS_DEFAULT)
      LibmeshPetscCallA(problem.comm().get(), KSPSetPCSide(ksp, getPetscPCSide(sys.getPCSide())));
  }
}
 
void
petscSetKSPDefaults(FEProblemBase & problem, KSP ksp)
{
  auto & es = problem.es();
 
  PetscReal rtol = es.parameters.get<Real>("linear solver tolerance");
  PetscReal atol = es.parameters.get<Real>("linear solver absolute tolerance");
 
  // MOOSE defaults this to -1 for some dumb reason
  if (atol < 0)
    atol = 1e-50;
 
  PetscReal maxits = es.parameters.get<unsigned int>("linear solver maximum iterations");
 
  // 1e100 is because we don't use divtol currently
  LibmeshPetscCallA(problem.comm().get(), KSPSetTolerances(ksp, rtol, atol, 1e100, maxits));
 
  petscSetDefaultPCSide(problem, ksp);
 
  petscSetDefaultKSPNormType(problem, ksp);
}
 
void
petscSetDefaults(FEProblemBase & problem)
{
  // Apply matrix-type options once the per-system matrix prefixes are known. This is different
  // from vectors: libMesh/PETSc vector construction already sees a global '-vec_type' option,
  // but prefixed matrix options such as '-nl0_mat_type' cannot match anything until we set the
  // matrix prefix here. Without this, a matrix may be constructed with the default type and keep
  // that type for the rest of the solve, unless we not only set the options prefix but also apply
  // the options to the matrix in this function call.
  applyMatrixTypeOptions(problem);
 
  for (const auto nl_index : make_range(problem.numNonlinearSystems()))
  {
    NonlinearSystemBase & nl = problem.getNonlinearSystemBase(nl_index);
 
    // dig out PETSc solver
    auto * const petsc_solver = cast_ptr<PetscNonlinearSolver<Number> *>(nl.nonlinearSolver());
 
    // Ensure we properly prefix SNES which in turn prefixes its KSP
    const char * snes_prefix = nullptr;
    std::string snes_prefix_str;
    if (nl.system().prefix_with_name())
    {
      snes_prefix_str = nl.system().prefix();
      snes_prefix = snes_prefix_str.c_str();
    }
    SNES snes = petsc_solver->snes(snes_prefix);
    KSP ksp;
    LibmeshPetscCallA(nl.comm().get(), SNESGetKSP(snes, &ksp));
    LibmeshPetscCallA(nl.comm().get(), SNESSetMaxLinearSolveFailures(snes, 1000000));
    LibmeshPetscCallA(nl.comm().get(), SNESSetCheckJacobianDomainError(snes, PETSC_TRUE));
    LibmeshPetscCallA(
        nl.comm().get(),
        SNESSetConvergenceTest(snes, petscNonlinearConverged, &problem, LIBMESH_PETSC_NULLPTR));
 
    petscSetKSPDefaults(problem, ksp);
  }
 
  for (auto sys_index : make_range(problem.numLinearSystems()))
  {
    // dig out PETSc solver
    LinearSystem & lin_sys = problem.getLinearSystem(sys_index);
    auto & lm_lin_sys = lin_sys.linearImplicitSystem();
    auto * const petsc_solver =
        dynamic_cast<PetscLinearSolver<Number> *>(lm_lin_sys.get_linear_solver());
    // Ensure we properly prefix KSP
    if (lm_lin_sys.prefix_with_name())
      petsc_solver->init(lm_lin_sys.prefix().c_str());
    else
      petsc_solver->init();
    // The KSP call here would initialize without a prefix if we hadn't "manually" performed
    // initialization above
    KSP ksp = petsc_solver->ksp();
 
    if (problem.hasLinearConvergenceObjects())
      LibmeshPetscCallA(
          lin_sys.comm().get(),
          KSPSetConvergenceTest(ksp, petscLinearConverged, &problem, LIBMESH_PETSC_NULLPTR));
 
    // We dont set the KSP defaults here because they seem to clash with the linear solve parameters
    // set in FEProblemBase::solveLinearSystem
  }
}
 
void
processSingletonMooseWrappedOptions(FEProblemBase & fe_problem, const InputParameters & params)
{
  setSolveTypeFromParams(fe_problem, params);
  setLineSearchFromParams(fe_problem, params);
  setMFFDTypeFromParams(fe_problem, params);
}
 
#define checkPrefix(prefix)                                                                        \
  mooseAssert(prefix[0] == '-',                                                                    \
              "Leading prefix character must be a '-'. Current prefix is '" << prefix << "'");     \
  mooseAssert((prefix.size() == 1) || (prefix.back() == '_'),                                      \
              "Terminating prefix character must be a '_'. Current prefix is '" << prefix << "'"); \
  mooseAssert(MooseUtils::isAllLowercase(prefix), "PETSc prefixes should be all lower-case")
 
void
storePetscOptions(FEProblemBase & fe_problem,
                  const std::string & prefix,
                  const ParallelParamObject & param_object)
{
  const auto & params = param_object.parameters();
  processSingletonMooseWrappedOptions(fe_problem, params);
 
  // The parameters contained in the Action
  const auto & petsc_options = params.get<MultiMooseEnum>("petsc_options");
  const auto & petsc_pair_options =
      params.get<MooseEnumItem, std::string>("petsc_options_iname", "petsc_options_value");
 
  // A reference to the PetscOptions object that contains the settings that will be used in the
  // solve
  auto & po = fe_problem.getPetscOptions();
 
  // First process the single petsc options/flags
  addPetscFlagsToPetscOptions(petsc_options, prefix, param_object, po);
 
  // Then process the option-value pairs
  addPetscPairsToPetscOptions(
      petsc_pair_options, fe_problem.mesh().dimension(), prefix, param_object, po);
}
 
void
setSolveTypeFromParams(FEProblemBase & fe_problem, const InputParameters & params)
{
  // Note: Options set in the Preconditioner block will override those set in the Executioner block
  if (params.isParamValid("solve_type") && !params.isParamValid("_use_eigen_value"))
  {
    // Extract the solve type
    const std::string & solve_type = params.get<MooseEnum>("solve_type");
    for (const auto i : make_range(fe_problem.numNonlinearSystems()))
      fe_problem.solverParams(i)._type = Moose::stringToEnum<Moose::SolveType>(solve_type);
  }
}
 
void
setLineSearchFromParams(FEProblemBase & fe_problem, const InputParameters & params)
{
  // Note: Options set in the Preconditioner block will override those set in the Executioner block
  if (params.isParamValid("line_search"))
  {
    const auto & line_search = params.get<MooseEnum>("line_search");
    for (const auto i : make_range(fe_problem.numNonlinearSystems()))
      if (fe_problem.solverParams(i)._line_search == Moose::LS_INVALID || line_search != "default")
      {
        Moose::LineSearchType enum_line_search =
            Moose::stringToEnum<Moose::LineSearchType>(line_search);
        fe_problem.solverParams(i)._line_search = enum_line_search;
        if (enum_line_search == LS_CONTACT || enum_line_search == LS_PROJECT)
        {
          NonlinearImplicitSystem * nl_system = dynamic_cast<NonlinearImplicitSystem *>(
              &fe_problem.getNonlinearSystemBase(i).system());
          if (!nl_system)
            mooseError("You've requested a line search but you must be solving an EigenProblem. "
                       "These two things are not consistent.");
          PetscNonlinearSolver<Real> * petsc_nonlinear_solver =
              dynamic_cast<PetscNonlinearSolver<Real> *>(nl_system->nonlinear_solver.get());
          if (!petsc_nonlinear_solver)
            mooseError("Currently the MOOSE line searches all use Petsc, so you "
                       "must use Petsc as your non-linear solver.");
          petsc_nonlinear_solver->linesearch_object =
              std::make_unique<ComputeLineSearchObjectWrapper>(fe_problem);
        }
      }
  }
}
 
void
setMFFDTypeFromParams(FEProblemBase & fe_problem, const InputParameters & params)
{
  if (params.isParamValid("mffd_type"))
  {
    const auto & mffd_type = params.get<MooseEnum>("mffd_type");
    for (const auto i : make_range(fe_problem.numNonlinearSystems()))
      fe_problem.solverParams(i)._mffd_type = Moose::stringToEnum<Moose::MffdType>(mffd_type);
  }
}
 
template <typename T>
void
checkUserProvidedPetscOption(const T & option, const ParallelParamObject & param_object)
{
  const auto & string_option = static_cast<const std::string &>(option);
  if (string_option[0] != '-')
    param_object.mooseError("PETSc option '", string_option, "' does not begin with '-'");
}
 
void
addPetscFlagsToPetscOptions(const MultiMooseEnum & petsc_flags,
                            std::string prefix,
                            const ParallelParamObject & param_object,
                            PetscOptions & po)
{
  prefix.insert(prefix.begin(), '-');
  checkPrefix(prefix);
 
  // Update the PETSc single flags
  for (const auto & option : petsc_flags)
  {
    checkUserProvidedPetscOption(option, param_object);
 
    const std::string & string_option = option.name();
 
    if (option == "-log_summary" || option == "-log_view")
      mooseError("The PETSc option \"-log_summary\" or \"-log_view\" can only be used on the "
                 "command line.  Please "
                 "remove it from the input file");
 
    // Update the stored items, but do not create duplicates
    const std::string prefixed_option = prefix + string_option.substr(1);
    if (!po.flags.isValueSet(prefixed_option))
    {
      po.flags.setAdditionalValue(prefixed_option);
      po.user_set_options.setAdditionalValue(prefixed_option);
    }
  }
}
 
void
setConvergedReasonFlags(FEProblemBase & fe_problem, std::string prefix)
{
  prefix.insert(prefix.begin(), '-');
  checkPrefix(prefix);
  libmesh_ignore(fe_problem); // avoid unused warnings for old PETSc
 
#if !PETSC_VERSION_LESS_THAN(3, 14, 0)
  // the boolean in these pairs denote whether the user has specified any of the reason flags in the
  // input file
  std::array<std::string, 2> reason_flags = {{"snes_converged_reason", "ksp_converged_reason"}};
 
  auto & po = fe_problem.getPetscOptions();
 
  for (const auto & reason_flag : reason_flags)
  {
    const auto full_flag = prefix + reason_flag;
    if (!po.flags.isValueSet(full_flag) && !po.dont_add_these_options.contains(full_flag) &&
        (std::find_if(po.pairs.begin(),
                      po.pairs.end(),
                      [&full_flag](auto & pair)
                      { return pair.first == (full_flag); }) == po.pairs.end()))
      po.pairs.emplace_back(full_flag, "::failed");
  }
#endif
}
 
void
addPetscPairsToPetscOptions(
    const std::vector<std::pair<MooseEnumItem, std::string>> & petsc_pair_options,
    const unsigned int mesh_dimension,
    std::string prefix,
    const ParallelParamObject & param_object,
    PetscOptions & po)
{
  prefix.insert(prefix.begin(), '-');
  checkPrefix(prefix);
 
  // Setup the name value pairs
  bool boomeramg_found = false;
  bool strong_threshold_found = false;
#if !PETSC_VERSION_LESS_THAN(3, 7, 0)
  bool superlu_dist_found = false;
  bool fact_pattern_found = false;
  bool tiny_pivot_found = false;
#endif
  std::string pc_description = "";
#if !PETSC_VERSION_LESS_THAN(3, 12, 0)
  // If users use HMG, we would like to set
  bool hmg_found = false;
  bool matptap_found = false;
  bool hmg_strong_threshold_found = false;
#endif
  std::vector<std::pair<std::string, std::string>> new_options;
 
  for (const auto & [option_name, option_value] : petsc_pair_options)
  {
    checkUserProvidedPetscOption(option_name, param_object);
 
    new_options.clear();
    const std::string prefixed_option_name =
        prefix + static_cast<const std::string &>(option_name).substr(1);
 
    // Do not add duplicate settings
    if (auto it =
            MooseUtils::findPair(po.pairs, po.pairs.begin(), prefixed_option_name, MooseUtils::Any);
        it == po.pairs.end())
    {
#if !PETSC_VERSION_LESS_THAN(3, 9, 0)
      if (option_name == "-pc_factor_mat_solver_package")
        new_options.emplace_back(prefix + "pc_factor_mat_solver_type", option_value);
#else
      if (option_name == "-pc_factor_mat_solver_type")
        new_options.push_back(prefix + "pc_factor_mat_solver_package", option_value);
#endif
 
      // Look for a pc description
      if (option_name == "-pc_type" || option_name == "-sub_pc_type" ||
          option_name == "-pc_hypre_type")
        pc_description += option_value + ' ';
 
#if !PETSC_VERSION_LESS_THAN(3, 12, 0)
      if (option_name == "-pc_type" && option_value == "hmg")
        hmg_found = true;
 
      // MPIAIJ for PETSc 3.12.0: -matptap_via
      // MAIJ for PETSc 3.12.0: -matmaijptap_via
      // MPIAIJ for PETSc 3.13 to 3.16: -matptap_via, -matproduct_ptap_via
      // MAIJ for PETSc 3.13 to 3.16: -matproduct_ptap_via
      // MPIAIJ for PETSc 3.17 and higher: -matptap_via, -mat_product_algorithm
      // MAIJ for PETSc 3.17 and higher: -mat_product_algorithm
#if !PETSC_VERSION_LESS_THAN(3, 17, 0)
      if (hmg_found && (option_name == "-matptap_via" || option_name == "-matmaijptap_via" ||
                        option_name == "-matproduct_ptap_via"))
        new_options.emplace_back(prefix + "mat_product_algorithm", option_value);
#elif !PETSC_VERSION_LESS_THAN(3, 13, 0)
      if (hmg_found && (option_name == "-matptap_via" || option_name == "-matmaijptap_via"))
        new_options.emplace_back(prefix + "matproduct_ptap_via", option_value);
#else
      if (hmg_found && (option_name == "-matproduct_ptap_via"))
      {
        new_options.emplace_back(prefix + "matptap_via", option_value);
        new_options.emplace_back(prefix + "matmaijptap_via", option_value);
      }
#endif
 
      if (option_name == "-matptap_via" || option_name == "-matmaijptap_via" ||
          option_name == "-matproduct_ptap_via" || option_name == "-mat_product_algorithm")
        matptap_found = true;
 
      // For 3D problems, we need to set this 0.7
      if (option_name == "-hmg_inner_pc_hypre_boomeramg_strong_threshold")
        hmg_strong_threshold_found = true;
#endif
      // This special case is common enough that we'd like to handle it for the user.
      if (option_name == "-pc_hypre_type" && option_value == "boomeramg")
        boomeramg_found = true;
      if (option_name == "-pc_hypre_boomeramg_strong_threshold")
        strong_threshold_found = true;
#if !PETSC_VERSION_LESS_THAN(3, 7, 0)
      if ((option_name == "-pc_factor_mat_solver_package" ||
           option_name == "-pc_factor_mat_solver_type") &&
          option_value == "superlu_dist")
        superlu_dist_found = true;
      if (option_name == "-mat_superlu_dist_fact")
        fact_pattern_found = true;
      if (option_name == "-mat_superlu_dist_replacetinypivot")
        tiny_pivot_found = true;
#endif
 
      if (!new_options.empty())
      {
        std::copy(new_options.begin(), new_options.end(), std::back_inserter(po.pairs));
        for (const auto & option : new_options)
          po.user_set_options.setAdditionalValue(option.first);
      }
      else
      {
        po.pairs.push_back(std::make_pair(prefixed_option_name, option_value));
        po.user_set_options.setAdditionalValue(prefixed_option_name);
      }
    }
    else
    {
      do
      {
        it->second = option_value;
        it = MooseUtils::findPair(po.pairs, std::next(it), prefixed_option_name, MooseUtils::Any);
      } while (it != po.pairs.end());
    }
  }
 
  // When running a 3D mesh with boomeramg, it is almost always best to supply a strong threshold
  // value. We will provide that for the user here if they haven't supplied it themselves.
  if (boomeramg_found && !strong_threshold_found && mesh_dimension == 3)
  {
    po.pairs.emplace_back(prefix + "pc_hypre_boomeramg_strong_threshold", "0.7");
    pc_description += "strong_threshold: 0.7 (auto)";
  }
 
#if !PETSC_VERSION_LESS_THAN(3, 12, 0)
  if (hmg_found && !hmg_strong_threshold_found && mesh_dimension == 3)
  {
    po.pairs.emplace_back(prefix + "hmg_inner_pc_hypre_boomeramg_strong_threshold", "0.7");
    pc_description += "strong_threshold: 0.7 (auto)";
  }
 
  // Default PETSc PtAP takes too much memory, and it is not quite useful
  // Let us switch to use new algorithm
  if (hmg_found && !matptap_found)
  {
#if !PETSC_VERSION_LESS_THAN(3, 17, 0)
    po.pairs.emplace_back(prefix + "mat_product_algorithm", "allatonce");
#elif !PETSC_VERSION_LESS_THAN(3, 13, 0)
    po.pairs.emplace_back(prefix + "matproduct_ptap_via", "allatonce");
#else
    po.pairs.emplace_back(prefix + "matptap_via", "allatonce");
    po.pairs.emplace_back(prefix + "matmaijptap_via", "allatonce");
#endif
  }
#endif
 
#if !PETSC_VERSION_LESS_THAN(3, 7, 0)
  // In PETSc-3.7.{0--4}, there is a bug when using superlu_dist, and we have to use
  // SamePattern_SameRowPerm, otherwise we use whatever we have in PETSc
  if (superlu_dist_found && !fact_pattern_found)
  {
    po.pairs.emplace_back(prefix + "mat_superlu_dist_fact",
#if PETSC_VERSION_LESS_THAN(3, 7, 5)
                          "SamePattern_SameRowPerm");
    pc_description += "mat_superlu_dist_fact: SamePattern_SameRowPerm ";
#else
                          "SamePattern");
    pc_description += "mat_superlu_dist_fact: SamePattern ";
#endif
  }
 
  // restore this superlu  option
  if (superlu_dist_found && !tiny_pivot_found)
  {
    po.pairs.emplace_back(prefix + "mat_superlu_dist_replacetinypivot", "1");
    pc_description += " mat_superlu_dist_replacetinypivot: true ";
  }
#endif
  // Set Preconditioner description
  if (!pc_description.empty() && prefix.size() > 1)
    po.pc_description += "[" + prefix.substr(1, prefix.size() - 2) + "]: ";
  po.pc_description += pc_description;
}
 
std::set<std::string>
getPetscValidLineSearches()
{
  return {"default", "shell", "none", "basic", "l2", "bt", "cp"};
}
 
InputParameters
getPetscValidParams()
{
  InputParameters params = emptyInputParameters();
 
  MooseEnum solve_type("PJFNK JFNK NEWTON FD LINEAR");
  params.addParam<MooseEnum>("solve_type",
                             solve_type,
                             "PJFNK: Preconditioned Jacobian-Free Newton Krylov "
                             "JFNK: Jacobian-Free Newton Krylov "
                             "NEWTON: Full Newton Solve "
                             "FD: Use finite differences to compute Jacobian "
                             "LINEAR: Solving a linear problem");
 
  MooseEnum mffd_type("wp ds", "wp");
  params.addParam<MooseEnum>("mffd_type",
                             mffd_type,
                             "Specifies the finite differencing type for "
                             "Jacobian-free solve types. Note that the "
                             "default is wp (for Walker and Pernice).");
 
  params.addParam<MultiMooseEnum>(
      "petsc_options", getCommonPetscFlags(), "Singleton PETSc options");
  params.addParam<MultiMooseEnum>(
      "petsc_options_iname", getCommonPetscKeys(), "Names of PETSc name/value pairs");
  params.addParam<std::vector<std::string>>(
      "petsc_options_value",
      "Values of PETSc name/value pairs (must correspond with \"petsc_options_iname\"");
  params.addParamNamesToGroup("solve_type petsc_options petsc_options_iname petsc_options_value "
                              "mffd_type",
                              "PETSc");
 
  return params;
}
 
MultiMooseEnum
getCommonSNESFlags()
{
  return MultiMooseEnum(
      "-ksp_monitor_snes_lg -snes_ksp_ew -snes_converged_reason "
      "-snes_ksp -snes_linesearch_monitor -snes_mf -snes_mf_operator -snes_monitor "
      "-snes_test_display -snes_view -snes_monitor_cancel",
      "",
      true);
}
 
MultiMooseEnum
getCommonKSPFlags()
{
  return MultiMooseEnum(
      "-ksp_converged_reason -ksp_gmres_modifiedgramschmidt -ksp_monitor", "", true);
}
 
MultiMooseEnum
getCommonPetscFlags()
{
  auto options = MultiMooseEnum("-dm_moose_print_embedding -dm_view", "", true);
  options.addValidName(getCommonKSPFlags());
  options.addValidName(getCommonSNESFlags());
  return options;
}
 
MultiMooseEnum
getCommonSNESKeys()
{
  return MultiMooseEnum("-snes_atol -snes_linesearch_type -snes_ls -snes_max_it -snes_rtol "
                        "-snes_divergence_tolerance -snes_type",
                        "",
                        true);
}
 
MultiMooseEnum
getCommonKSPKeys()
{
  return MultiMooseEnum("-ksp_atol -ksp_gmres_restart -ksp_max_it -ksp_pc_side -ksp_rtol "
                        "-ksp_type -sub_ksp_type",
                        "",
                        true);
}
MultiMooseEnum
getCommonPetscKeys()
{
  auto options = MultiMooseEnum("-mat_fd_coloring_err -mat_fd_type -mat_mffd_type "
                                "-pc_asm_overlap -pc_factor_levels "
                                "-pc_factor_mat_ordering_type -pc_hypre_boomeramg_grid_sweeps_all "
                                "-pc_hypre_boomeramg_max_iter "
                                "-pc_hypre_boomeramg_strong_threshold -pc_hypre_type -pc_type "
                                "-sub_pc_type",
                                "",
                                true);
  options.addValidName(getCommonKSPKeys());
  options.addValidName(getCommonSNESKeys());
  return options;
}
 
bool
isSNESVI(FEProblemBase & fe_problem)
{
  const PetscOptions & petsc = fe_problem.getPetscOptions();
 
  int argc;
  char ** args;
  LibmeshPetscCallA(fe_problem.comm().get(), PetscGetArgs(&argc, &args));
 
  std::vector<std::string> cml_arg;
  for (int i = 0; i < argc; i++)
    cml_arg.push_back(args[i]);
 
  if (MooseUtils::findPair(petsc.pairs, petsc.pairs.begin(), MooseUtils::Any, "vinewtonssls") ==
          petsc.pairs.end() &&
      MooseUtils::findPair(petsc.pairs, petsc.pairs.begin(), MooseUtils::Any, "vinewtonrsls") ==
          petsc.pairs.end() &&
      std::find(cml_arg.begin(), cml_arg.end(), "vinewtonssls") == cml_arg.end() &&
      std::find(cml_arg.begin(), cml_arg.end(), "vinewtonrsls") == cml_arg.end())
    return false;
 
  return true;
}
 
void
setSinglePetscOption(const std::string & name,
                     const std::string & value /*=""*/,
                     FEProblemBase * const problem /*=nullptr*/)
{
  static const TIMPI::Communicator comm_world(PETSC_COMM_WORLD);
  const TIMPI::Communicator & comm = problem ? problem->comm() : comm_world;
  LibmeshPetscCallA(comm.get(),
                    PetscOptionsSetValue(LIBMESH_PETSC_NULLPTR,
                                         name.c_str(),
                                         value == "" ? LIBMESH_PETSC_NULLPTR : value.c_str()));
  // Create a single option data base so that we can use PETSC's internal option checking which
  // is case insensitive. This is better than re-implementing case-insensitive checks here in this
  // TU
  ::PetscOptions single_option;
  LibmeshPetscCallA(comm.get(), PetscOptionsCreate(&single_option));
  LibmeshPetscCallA(comm.get(),
                    PetscOptionsSetValue(single_option,
                                         name.c_str(),
                                         value == "" ? LIBMESH_PETSC_NULLPTR : value.c_str()));
  auto check_problem = [problem, &name]()
  {
    if (!problem)
      mooseError(
          "Setting the option '",
          name,
          "' requires passing a 'problem' parameter. Contact a developer of your application "
          "to have them update their code. If in doubt, reach out to the MOOSE team on Github "
          "discussions");
  };
 
  // Select vector type from user-passed PETSc options
  if (petscOptionsHasName(single_option, "-vec_type"))
  {
    check_problem();
    applyVectorTypeOptions(*problem);
  }
  // First do a cheap suffix check so unrelated PETSc options do not pay for looping over every
  // solver system. Once we know the name looks like a matrix-type option, rely on PETSc's
  // option lookup for the actual case-insensitive and prefix-aware matching.
  else if (problem && mightBeMatTypeOption(name))
  {
    errorOnUnprefixedMatTypeOption(single_option, *problem);
 
    for (const auto i : index_range(problem->_solver_systems))
    {
      const auto & solver_sys_name = problem->_solver_sys_names[i];
      if (!petscOptionsHasName(single_option, "-mat_type", solver_sys_name + "_"))
        continue;
 
      if (problem->solverParams(i)._type == Moose::ST_JFNK)
        mooseError("Setting option '", name, "' is incompatible with a JFNK 'solve_type'");
 
      applyMatrixTypeOptions(*problem, i);
      break;
    }
  }
 
  LibmeshPetscCallA(comm.get(), PetscOptionsDestroy(&single_option));
}
 
void
setSinglePetscOptionIfAppropriate(const MultiMooseEnum & dont_add_these_options,
                                  const std::string & name,
                                  const std::string & value /*=""*/,
                                  FEProblemBase * const problem /*=nullptr*/)
{
  if (!dont_add_these_options.contains(name))
    setSinglePetscOption(name, value, problem);
}
 
void
colorAdjacencyMatrix(PetscScalar * adjacency_matrix,
                     unsigned int size,
                     unsigned int colors,
                     std::vector<unsigned int> & vertex_colors,
                     const char * coloring_algorithm)
{
  // Mat A will be a dense matrix from the incoming data structure
  Mat A;
  LibmeshPetscCallA(PETSC_COMM_SELF, MatCreate(PETSC_COMM_SELF, &A));
  LibmeshPetscCallA(PETSC_COMM_SELF, MatSetSizes(A, size, size, size, size));
  LibmeshPetscCallA(PETSC_COMM_SELF, MatSetType(A, MATSEQDENSE));
  // PETSc requires a non-const data array to populate the matrix
  LibmeshPetscCallA(PETSC_COMM_SELF, MatSeqDenseSetPreallocation(A, adjacency_matrix));
  LibmeshPetscCallA(PETSC_COMM_SELF, MatAssemblyBegin(A, MAT_FINAL_ASSEMBLY));
  LibmeshPetscCallA(PETSC_COMM_SELF, MatAssemblyEnd(A, MAT_FINAL_ASSEMBLY));
 
  // Convert A to a sparse matrix
#if PETSC_VERSION_LESS_THAN(3, 7, 0)
  LibmeshPetscCallA(PETSC_COMM_SELF, MatConvert(A, MATAIJ, MAT_REUSE_MATRIX, &A));
#else
  LibmeshPetscCallA(PETSC_COMM_SELF, MatConvert(A, MATAIJ, MAT_INPLACE_MATRIX, &A));
#endif
 
  ISColoring iscoloring;
  MatColoring mc;
  LibmeshPetscCallA(PETSC_COMM_SELF, MatColoringCreate(A, &mc));
  LibmeshPetscCallA(PETSC_COMM_SELF, MatColoringSetType(mc, coloring_algorithm));
  LibmeshPetscCallA(PETSC_COMM_SELF, MatColoringSetMaxColors(mc, static_cast<PetscInt>(colors)));
 
  // Petsc normally colors by distance two (neighbors of neighbors), we just want one
  LibmeshPetscCallA(PETSC_COMM_SELF, MatColoringSetDistance(mc, 1));
  LibmeshPetscCallA(PETSC_COMM_SELF, MatColoringSetFromOptions(mc));
  LibmeshPetscCallA(PETSC_COMM_SELF, MatColoringApply(mc, &iscoloring));
 
  PetscInt nn;
  IS * is;
#if PETSC_RELEASE_LESS_THAN(3, 12, 0)
  LibmeshPetscCallA(PETSC_COMM_SELF, ISColoringGetIS(iscoloring, &nn, &is));
#else
  LibmeshPetscCallA(PETSC_COMM_SELF, ISColoringGetIS(iscoloring, PETSC_USE_POINTER, &nn, &is));
#endif
 
  if (nn > static_cast<PetscInt>(colors))
    throw std::runtime_error("Not able to color with designated number of colors");
 
  for (int i = 0; i < nn; i++)
  {
    PetscInt isize;
    const PetscInt * indices;
    LibmeshPetscCallA(PETSC_COMM_SELF, ISGetLocalSize(is[i], &isize));
    LibmeshPetscCallA(PETSC_COMM_SELF, ISGetIndices(is[i], &indices));
    for (int j = 0; j < isize; j++)
    {
      mooseAssert(indices[j] < static_cast<PetscInt>(vertex_colors.size()), "Index out of bounds");
      vertex_colors[indices[j]] = i;
    }
    LibmeshPetscCallA(PETSC_COMM_SELF, ISRestoreIndices(is[i], &indices));
  }
 
  LibmeshPetscCallA(PETSC_COMM_SELF, MatDestroy(&A));
  LibmeshPetscCallA(PETSC_COMM_SELF, MatColoringDestroy(&mc));
  LibmeshPetscCallA(PETSC_COMM_SELF, ISColoringDestroy(&iscoloring));
}
 
void
dontAddPetscFlag(const std::string & flag, PetscOptions & petsc_options)
{
  if (!petsc_options.dont_add_these_options.contains(flag))
    petsc_options.dont_add_these_options.setAdditionalValue(flag);
}
 
void
dontAddNonlinearConvergedReason(FEProblemBase & fe_problem)
{
  dontAddPetscFlag("-snes_converged_reason", fe_problem.getPetscOptions());
}
 
void
dontAddLinearConvergedReason(FEProblemBase & fe_problem)
{
  dontAddPetscFlag("-ksp_converged_reason", fe_problem.getPetscOptions());
}
 
void
dontAddCommonKSPOptions(FEProblemBase & fe_problem)
{
  auto & petsc_options = fe_problem.getPetscOptions();
  for (const auto & flag : getCommonKSPFlags().getNames())
    dontAddPetscFlag(flag, petsc_options);
  for (const auto & key : getCommonKSPKeys().getNames())
    dontAddPetscFlag(key, petsc_options);
}
 
void
dontAddCommonSNESOptions(FEProblemBase & fe_problem)
{
  dontAddCommonSNESOptions(fe_problem, "");
}
 
void
dontAddCommonSNESOptions(FEProblemBase & fe_problem, const std::string & prefix)
{
  auto & petsc_options = fe_problem.getPetscOptions();
  for (const auto & flag : getCommonSNESFlags().getNames())
    dontAddPetscFlag("-" + prefix + flag.substr(1), petsc_options);
  for (const auto & key : getCommonSNESKeys().getNames())
    dontAddPetscFlag("-" + prefix + key.substr(1), petsc_options);
}
 
std::unique_ptr<PetscMatrix<Number>>
createMatrixFromFile(const libMesh::Parallel::Communicator & comm,
                     Mat & mat,
                     const std::string & binary_mat_file,
                     const unsigned int mat_number_to_load)
{
  LibmeshPetscCallA(comm.get(), MatCreate(comm.get(), &mat));
  PetscViewer matviewer;
  LibmeshPetscCallA(
      comm.get(),
      PetscViewerBinaryOpen(comm.get(), binary_mat_file.c_str(), FILE_MODE_READ, &matviewer));
  for (unsigned int i = 0; i < mat_number_to_load; ++i)
    LibmeshPetscCallA(comm.get(), MatLoad(mat, matviewer));
  LibmeshPetscCallA(comm.get(), PetscViewerDestroy(&matviewer));
 
  return std::make_unique<PetscMatrix<Number>>(mat, comm);
}
 
void
registerPetscCitation(const std::string & bibtex)
{
  // PETSc concatenates registered entries verbatim, so ensure a trailing newline to keep entries
  // separated. A null "set" flag registers unconditionally; callers deduplicate by citation key.
  const std::string entry = (!bibtex.empty() && bibtex.back() == '\n') ? bibtex : bibtex + "\n";
  LibmeshPetscCallA(PETSC_COMM_WORLD, PetscCitationsRegister(entry.c_str(), nullptr));
}
 
} // Namespace PetscSupport
} // Namespace MOOSE