Moose::SlepcSupport Namespace Reference

Functions
InputParameters	getSlepcValidParams (InputParameters &params)
InputParameters	getSlepcEigenProblemValidParams ()
	Retrieve valid params that allow users to specify eigen problem configuration.
void	storeSolveType (FEProblemBase &fe_problem, const InputParameters &params)
	Set solve type into eigen problem (solverParams)
void	setEigenProblemSolverParams (EigenProblem &eigen_problem, const InputParameters &params)
	Retrieve eigen problem params from 'params', and then set these params into SolverParams.
void	slepcSetOptions (EigenProblem &eigen_problem, SolverParams &solver_params, const InputParameters &params)
	Push all SLEPc/PETSc options into SLEPc/PETSc side.
void	setSlepcEigenSolverTolerances (EigenProblem &eigen_problem, const SolverParams &solver_params, const InputParameters &params)
	Control eigen solver tolerances via SLEPc options.
void	setFreeNonlinearPowerIterations (unsigned int free_power_iterations)
	Set SLEPc/PETSc options to trigger free power iteration.
void	clearFreeNonlinearPowerIterations (const InputParameters &params)
void	moosePetscSNESFormMatrixTag (SNES snes, Vec x, Mat eigen_mat, SparseMatrix< Number > &all_dofs_mat, void *ctx, TagID tag)
	Form matrix according to tag.
void	moosePetscSNESFormMatricesTags (SNES snes, Vec x, std::vector< Mat > &eigen_mats, std::vector< SparseMatrix< Number > * > &all_dofs_mats, void *ctx, const std::set< TagID > &tags)
	Form multiple matrices for multiple tags.
PetscErrorCode	mooseSlepcEigenFormJacobianA (SNES snes, Vec x, Mat jac, Mat pc, void *ctx)
	Form Jacobian matrix A.
PetscErrorCode	mooseSlepcEigenFormJacobianB (SNES snes, Vec x, Mat jac, Mat pc, void *ctx)
	Form Jacobian matrix B.
PetscErrorCode	mooseSlepcEigenFormFunctionA (SNES snes, Vec x, Vec r, void *ctx)
	Form function residual Ax.
PetscErrorCode	mooseSlepcEigenFormFunctionB (SNES snes, Vec x, Vec r, void *ctx)
	Form function residual Bx.
PetscErrorCode	mooseSlepcEigenFormFunctionAB (SNES snes, Vec x, Vec Ax, Vec Bx, void *ctx)
	Form function residual Ax-Bx.
PetscErrorCode	mooseSlepcStoppingTest (EPS eps, PetscInt its, PetscInt max_it, PetscInt nconv, PetscInt nev, EPSConvergedReason *reason, void *ctx)
	A customized convergence checker.
PetscErrorCode	mooseSlepcEPSGetSNES (EPS eps, SNES *snes)
	Retrieve SNES from EPS.
PetscErrorCode	mooseSlepcEPSMonitor (EPS eps, PetscInt its, PetscInt nconv, PetscScalar *eigr, PetscScalar *eigi, PetscReal *errest, PetscInt nest, void *mctx)
	A customized solver monitor to print out eigenvalue.
PetscErrorCode	mooseSlepcEPSSNESSetUpOptionPrefix (EPS eps)
	Get rid of prefix "-eps_power" for SNES, KSP, PC, etc.
PetscErrorCode	mooseSlepcEPSSNESSetCustomizePC (EPS eps)
	Attach a customized PC.
PetscErrorCode	mooseSlepcEPSSNESKSPSetPCSide (FEProblemBase &problem, EPS eps)
	Allow users to specify PC side.
void	attachCallbacksToMat (EigenProblem &eigen_problem, Mat mat, bool eigen)
	Attach call backs to mat.
PetscErrorCode	mooseMatMult_Eigen (Mat mat, Vec x, Vec y)
	Implement MatMult via function evaluation for Bx.
PetscErrorCode	mooseMatMult_NonEigen (Mat mat, Vec x, Vec y)
	Implement MatMult via function evaluation for Ax.
void	setOperationsForShellMat (EigenProblem &eigen_problem, Mat mat, bool eigen)
	Set operations to shell mat.
PETSC_EXTERN PetscErrorCode	PCCreate_MoosePC (PC pc)
	Create a preconditioner from moose side.
PETSC_EXTERN PetscErrorCode	registerPCToPETSc ()
	Let PETSc know there is a preconditioner.
PetscErrorCode	PCDestroy_MoosePC (PC pc)
	Destroy preconditioner.
PetscErrorCode	PCView_MoosePC (PC pc, PetscViewer viewer)
	View preconditioner.
PetscErrorCode	PCApply_MoosePC (PC pc, Vec x, Vec y)
	Preconditioner application.
PetscErrorCode	PCSetUp_MoosePC (PC pc)
	Setup preconditioner.
PetscErrorCode	mooseSlepcEigenFormFunctionMFFD (void *ctx, Vec x, Vec r)
	Function call for MFFD.
void	setEigenProblemOptions (SolverParams &solver_params, const MultiMooseEnum &dont_add_these_options)
void	setWhichEigenPairsOptions (SolverParams &solver_params, const MultiMooseEnum &dont_add_these_options)
void	setNewtonPetscOptions (SolverParams &solver_params, const InputParameters &params)
void	setNonlinearPowerOptions (SolverParams &solver_params)
void	setEigenSolverOptions (SolverParams &solver_params, const InputParameters &params)
PetscErrorCode	mooseEPSFormMatrices (EigenProblem &eigen_problem, EPS eps, Vec x, void *ctx)
void	moosePetscSNESFormFunction (SNES, Vec x, Vec r, void *ctx, TagID tag)
PetscErrorCode	mooseSlepcEigenFormNorm (SNES, Vec, PetscReal *norm, void *ctx)

Variables
const int	subspace_factor = 2

Function Documentation

attachCallbacksToMat()

void Moose::SlepcSupport::attachCallbacksToMat	(	EigenProblem &	eigen_problem,
		Mat	mat,
		bool	eigen
	)

Attach call backs to mat.

SLEPc solver will retrieve callbacks we attached here

Definition at line 1060 of file SlepcSupport.C.

{
  // Recall that we are solving the potentially nonlinear problem:
  // F(x) = A(x) - \lambda B(x) = 0
  //
  // To solve this, we can use Newton's method: J \Delta x = -F
  // Generally we will approximate J using matrix free methods. However, in order to solve the
  // linearized system efficiently, we typically will need preconditioning. Typically we will build
  // the preconditioner only from A, but we also have the option to include information from B
 
  // Attach the Jacobian computation function. If \p mat is the "eigen" matrix corresponding to B,
  // then attach our JacobianB computation routine, else the matrix corresponds to A, and we attach
  // the JacobianA computation routine
  LibmeshPetscCallA(
      eigen_problem.comm().get(),
      PetscObjectComposeFunction((PetscObject)mat,
                                 "formJacobian",
                                 eigen ? Moose::SlepcSupport::mooseSlepcEigenFormJacobianB
                                       : Moose::SlepcSupport::mooseSlepcEigenFormJacobianA));
 
  // Attach the residual computation function. If \p mat is the "eigen" matrix corresponding to B,
  // then attach our FunctionB computation routine, else the matrix corresponds to A, and we attach
  // the FunctionA computation routine
  LibmeshPetscCallA(
      eigen_problem.comm().get(),
      PetscObjectComposeFunction((PetscObject)mat,
                                 "formFunction",
                                 eigen ? Moose::SlepcSupport::mooseSlepcEigenFormFunctionB
                                       : Moose::SlepcSupport::mooseSlepcEigenFormFunctionA));
 
  // It's also beneficial to be able to evaluate both A and B residuals at once
  LibmeshPetscCallA(eigen_problem.comm().get(),
                    PetscObjectComposeFunction((PetscObject)mat,
                                               "formFunctionAB",
                                               Moose::SlepcSupport::mooseSlepcEigenFormFunctionAB));
 
  // Users may choose to provide a custom measure of the norm of B (Bx for a linear system)
  if (eigen_problem.bxNormProvided())
    LibmeshPetscCallA(eigen_problem.comm().get(),
                      PetscObjectComposeFunction((PetscObject)mat,
                                                 "formNorm",
                                                 Moose::SlepcSupport::mooseSlepcEigenFormNorm));
 
  // Finally we need to attach the "context" object, which is our EigenProblem, to the matrices so
  // that eventually when we get callbacks from SLEPc we can call methods on the EigenProblem
  PetscContainer container;
  LibmeshPetscCallA(eigen_problem.comm().get(),
                    PetscContainerCreate(eigen_problem.comm().get(), &container));
  LibmeshPetscCallA(eigen_problem.comm().get(),
                    PetscContainerSetPointer(container, &eigen_problem));
  LibmeshPetscCallA(
      eigen_problem.comm().get(),
      PetscObjectCompose((PetscObject)mat, "formJacobianCtx", (PetscObject)container));
  LibmeshPetscCallA(
      eigen_problem.comm().get(),
      PetscObjectCompose((PetscObject)mat, "formFunctionCtx", (PetscObject)container));
  if (eigen_problem.bxNormProvided())
    LibmeshPetscCallA(eigen_problem.comm().get(),
                      PetscObjectCompose((PetscObject)mat, "formNormCtx", (PetscObject)container));
 
  LibmeshPetscCallA(eigen_problem.comm().get(), PetscContainerDestroy(&container));
}

Referenced by NonlinearEigenSystem::attachSLEPcCallbacks().

clearFreeNonlinearPowerIterations()

void Moose::SlepcSupport::clearFreeNonlinearPowerIterations

(

const InputParameters &

params

)

Definition at line 417 of file SlepcSupport.C.

{
  Moose::PetscSupport::setSinglePetscOption("-eps_power_update", "1");
  Moose::PetscSupport::setSinglePetscOption("-eps_max_it", "1");
  Moose::PetscSupport::setSinglePetscOption("-snes_max_it",
                                            stringify(params.get<unsigned int>("nl_max_its")));
  Moose::PetscSupport::setSinglePetscOption("-snes_rtol",
                                            stringify(params.get<Real>("nl_rel_tol")));
  Moose::PetscSupport::setSinglePetscOption("-eps_tol", stringify(params.get<Real>("eigen_tol")));
}

Referenced by EigenProblem::doFreeNonlinearPowerIterations().

getSlepcEigenProblemValidParams()

InputParameters Moose::SlepcSupport::getSlepcEigenProblemValidParams

(

)

Retrieve valid params that allow users to specify eigen problem configuration.

Definition at line 64 of file SlepcSupport.C.

{
  InputParameters params = emptyInputParameters();
 
  // We are solving a Non-Hermitian eigenvalue problem by default
  MooseEnum eigen_problem_type("HERMITIAN NON_HERMITIAN GEN_HERMITIAN GEN_NON_HERMITIAN "
                               "GEN_INDEFINITE POS_GEN_NON_HERMITIAN SLEPC_DEFAULT",
                               "GEN_NON_HERMITIAN");
  params.addParam<MooseEnum>(
      "eigen_problem_type",
      eigen_problem_type,
      "Type of the eigenvalue problem we are solving "
      "HERMITIAN: Hermitian "
      "NON_HERMITIAN: Non-Hermitian "
      "GEN_HERMITIAN: Generalized Hermitian "
      "GEN_NON_HERMITIAN: Generalized Non-Hermitian "
      "GEN_INDEFINITE: Generalized indefinite Hermitian "
      "POS_GEN_NON_HERMITIAN: Generalized Non-Hermitian with positive (semi-)definite B "
      "SLEPC_DEFAULT: Use whatever SLEPC has by default ");
 
  // Which eigenvalues are we interested in
  MooseEnum which_eigen_pairs("LARGEST_MAGNITUDE SMALLEST_MAGNITUDE LARGEST_REAL SMALLEST_REAL "
                              "LARGEST_IMAGINARY SMALLEST_IMAGINARY TARGET_MAGNITUDE TARGET_REAL "
                              "TARGET_IMAGINARY ALL_EIGENVALUES SLEPC_DEFAULT");
  params.addParam<MooseEnum>("which_eigen_pairs",
                             which_eigen_pairs,
                             "Which eigenvalue pairs to obtain from the solution "
                             "LARGEST_MAGNITUDE "
                             "SMALLEST_MAGNITUDE "
                             "LARGEST_REAL "
                             "SMALLEST_REAL "
                             "LARGEST_IMAGINARY "
                             "SMALLEST_IMAGINARY "
                             "TARGET_MAGNITUDE "
                             "TARGET_REAL "
                             "TARGET_IMAGINARY "
                             "ALL_EIGENVALUES "
                             "SLEPC_DEFAULT ");
 
  params.addParam<unsigned int>("n_eigen_pairs", 1, "The number of eigen pairs");
  params.addParam<unsigned int>("n_basis_vectors", 3, "The dimension of eigen subspaces");
 
  params.addParam<Real>("eigen_tol", 1.0e-4, "Relative Tolerance for Eigen Solver");
  params.addParam<unsigned int>("eigen_max_its", 10000, "Max Iterations for Eigen Solver");
 
  params.addParam<Real>("l_abs_tol", 1e-50, "Absolute Tolerances for Linear Solver");
 
  params.addParam<unsigned int>("free_power_iterations", 4, "The number of free power iterations");
 
  params.addParam<unsigned int>(
      "extra_power_iterations", 0, "The number of extra free power iterations");
 
  params.addParamNamesToGroup(
      "eigen_problem_type which_eigen_pairs n_eigen_pairs n_basis_vectors eigen_tol eigen_max_its "
      "free_power_iterations extra_power_iterations",
      "Eigen Solver");
  params.addParamNamesToGroup("l_abs_tol", "Linear solver");
 
  return params;
}

Referenced by EigenProblemSolve::validParams().

getSlepcValidParams()

InputParameters Moose::SlepcSupport::getSlepcValidParams

(

InputParameters &

params

)

Returns

InputParameters object containing the SLEPC related parameters

The output of this function should be added to the the parameters object of the overarching class

See also

EigenProblem

Definition at line 35 of file SlepcSupport.C.

{
  MooseEnum solve_type("POWER ARNOLDI KRYLOVSCHUR JACOBI_DAVIDSON "
                       "NONLINEAR_POWER NEWTON PJFNK PJFNKMO JFNK",
                       "PJFNK");
  params.set<MooseEnum>("solve_type") = solve_type;
 
  params.setDocString("solve_type",
                      "POWER: Power / Inverse / RQI "
                      "ARNOLDI: Arnoldi "
                      "KRYLOVSCHUR: Krylov-Schur "
                      "JACOBI_DAVIDSON: Jacobi-Davidson "
                      "NONLINEAR_POWER: Nonlinear Power "
                      "NEWTON: Newton "
                      "PJFNK: Preconditioned Jacobian-free Newton-Kyrlov"
                      "JFNK: Jacobian-free Newton-Kyrlov"
                      "PJFNKMO: Preconditioned Jacobian-free Newton-Kyrlov with Matrix Only");
 
  // When the eigenvalue problems is reformed as a coupled nonlinear system,
  // we use part of Jacobian as the preconditioning matrix.
  // Because the difference between the Jacobian and the preconditioning matrix is not small,
  // the linear solver KSP can not reduce the residual much. After several tests,
  // we find 1e-2 is a reasonable choice.
  params.set<Real>("l_tol") = 1e-2;
 
  return params;
}

Referenced by EigenProblemSolve::validParams().

mooseEPSFormMatrices()

PetscErrorCode Moose::SlepcSupport::mooseEPSFormMatrices	(	EigenProblem &	eigen_problem,
		EPS	eps,
		Vec	x,
		void *	ctx
	)

Definition at line 571 of file SlepcSupport.C.

{
  ST st;
  Mat A, B;
  PetscBool aisshell, bisshell;
  PetscFunctionBegin;
 
  if (eigen_problem.constantMatrices() && eigen_problem.wereMatricesFormed())
    PetscFunctionReturn(PETSC_SUCCESS);
 
  if (eigen_problem.onLinearSolver())
    // We reach here during linear iteration when solve type is PJFNKMO.
    // We will use the matrices assembled at the beginning of this Newton
    // iteration for the following residual evaluation.
    PetscFunctionReturn(PETSC_SUCCESS);
 
  NonlinearEigenSystem & eigen_nl = eigen_problem.getCurrentNonlinearEigenSystem();
  auto & sys = eigen_nl.sys();
  SNES snes = eigen_nl.getSNES();
  // Rest ST state so that we can retrieve matrices
  LibmeshPetscCallQ(EPSGetST(eps, &st));
  LibmeshPetscCallQ(STResetMatrixState(st));
  LibmeshPetscCallQ(EPSGetOperators(eps, &A, &B));
  LibmeshPetscCallQ(PetscObjectTypeCompare((PetscObject)A, MATSHELL, &aisshell));
  LibmeshPetscCallQ(PetscObjectTypeCompare((PetscObject)B, MATSHELL, &bisshell));
  if (aisshell || bisshell)
  {
    SETERRQ(PetscObjectComm((PetscObject)eps),
            PETSC_ERR_ARG_INCOMP,
            "A and B matrices can not be shell matrices when using PJFNKMO \n");
  }
  // Form A and B
  std::vector<Mat> mats = {A, B};
  std::vector<SparseMatrix<Number> *> libmesh_mats = {&sys.get_matrix_A(), &sys.get_matrix_B()};
  moosePetscSNESFormMatricesTags(
      snes, x, mats, libmesh_mats, ctx, {eigen_nl.nonEigenMatrixTag(), eigen_nl.eigenMatrixTag()});
  eigen_problem.wereMatricesFormed(true);
  PetscFunctionReturn(PETSC_SUCCESS);
}

Referenced by mooseSlepcEigenFormFunctionA(), mooseSlepcEigenFormFunctionAB(), mooseSlepcEigenFormFunctionB(), and mooseSlepcEigenFormJacobianA().

mooseMatMult_Eigen()

PetscErrorCode Moose::SlepcSupport::mooseMatMult_Eigen	(	Mat	mat,
		Vec	x,
		Vec	y
	)

Implement MatMult via function evaluation for Bx.

Definition at line 1124 of file SlepcSupport.C.

{
  PetscFunctionBegin;
  void * ctx = nullptr;
  LibmeshPetscCallQ(MatShellGetContext(mat, &ctx));
 
  if (!ctx)
    mooseError("No context is set for shell matrix ");
 
  EigenProblem * eigen_problem = static_cast<EigenProblem *>(ctx);
  NonlinearEigenSystem & eigen_nl = eigen_problem->getCurrentNonlinearEigenSystem();
 
  evaluateResidual(*eigen_problem, x, r, eigen_nl.eigenVectorTag());
 
  if (eigen_problem->negativeSignEigenKernel())
    LibmeshPetscCallQ(VecScale(r, -1.));
 
  PetscFunctionReturn(PETSC_SUCCESS);
}

Referenced by setOperationsForShellMat().

mooseMatMult_NonEigen()

PetscErrorCode Moose::SlepcSupport::mooseMatMult_NonEigen	(	Mat	mat,
		Vec	x,
		Vec	y
	)

Implement MatMult via function evaluation for Ax.

Definition at line 1145 of file SlepcSupport.C.

{
  PetscFunctionBegin;
  void * ctx = nullptr;
  LibmeshPetscCallQ(MatShellGetContext(mat, &ctx));
 
  if (!ctx)
    mooseError("No context is set for shell matrix ");
 
  EigenProblem * eigen_problem = static_cast<EigenProblem *>(ctx);
  NonlinearEigenSystem & eigen_nl = eigen_problem->getCurrentNonlinearEigenSystem();
 
  evaluateResidual(*eigen_problem, x, r, eigen_nl.nonEigenVectorTag());
 
  PetscFunctionReturn(PETSC_SUCCESS);
}

Referenced by setOperationsForShellMat().

moosePetscSNESFormFunction()

void Moose::SlepcSupport::moosePetscSNESFormFunction	(	SNES	,
		Vec	x,
		Vec	r,
		void *	ctx,
		TagID	tag
	)

Definition at line 902 of file SlepcSupport.C.

{
  EigenProblem * eigen_problem = static_cast<EigenProblem *>(ctx);
  evaluateResidual(*eigen_problem, x, r, tag);
}

Referenced by mooseSlepcEigenFormFunctionA(), and mooseSlepcEigenFormFunctionB().

moosePetscSNESFormMatricesTags()

void Moose::SlepcSupport::moosePetscSNESFormMatricesTags	(	SNES	snes,
		Vec	x,
		std::vector< Mat > &	eigen_mats,
		std::vector< SparseMatrix< Number > * > &	all_dofs_mats,
		void *	ctx,
		const std::set< TagID > &	tags
	)

Form multiple matrices for multiple tags.

This function is conceptually the same as moosePetscSNESForMatrixTag() but allows passing in multiple matrices to be formed

Definition at line 710 of file SlepcSupport.C.

{
  EigenProblem * eigen_problem = static_cast<EigenProblem *>(ctx);
  auto & nl = eigen_problem->getCurrentNonlinearEigenSystem();
  auto & sys = nl.sys();
  auto & dof_map = sys.get_dof_map();
 
#ifndef NDEBUG
  for (const auto i : index_range(eigen_mats))
    mooseAssert(!dof_map.n_constrained_dofs() ==
                    (eigen_mats[i] == cast_ptr<PetscMatrix<Number> *>(all_dofs_mats[i])->mat()),
                "If we do not have constrained dofs, then mat and libmesh_mat should be the same. "
                "Conversely, if we do have constrained dofs, they must be different");
#endif
 
  updateCurrentLocalSolution(sys, x);
 
  for (auto * const all_dofs_mat : all_dofs_mats)
    if (!eigen_problem->constJacobian())
      all_dofs_mat->zero();
 
  eigen_problem->computeMatricesTags(*sys.current_local_solution.get(), all_dofs_mats, tags);
 
  if (dof_map.n_constrained_dofs())
    for (const auto i : index_range(eigen_mats))
    {
      PetscMatrix<Number> wrapped_eigen_mat(eigen_mats[i], sys.comm());
      sys.copy_super_to_sub(*all_dofs_mats[i], wrapped_eigen_mat);
    }
}

Referenced by mooseEPSFormMatrices(), and mooseSlepcEigenFormJacobianA().

moosePetscSNESFormMatrixTag()

void Moose::SlepcSupport::moosePetscSNESFormMatrixTag	(	SNES	snes,
		Vec	x,
		Mat	eigen_mat,
		SparseMatrix< Number > &	all_dofs_mat,
		void *	ctx,
		TagID	tag
	)

Form matrix according to tag.

Parameters

snes	The PETSc nonlinear solver context associated with the eigensolver
x	The current eigen solution vector from SLEPc. This vector will have a size corresponding to the number of unconstrained degrees of freedom
eigen_mat	The matrix we want to form coming from SLEPc. As for x, this will be sized corresponding to the number of unconstrained degrees of freedom
all_dofs_mat	This matrix corresponds to the "global" matrix, representing both unconstrained and constrained degrees of freedom. We will pass this matrix to MOOSE's assembly framework and then if the problem has constraints we will selectively copy the submatrix representing unconstrained dofs to eigen_mat. (Else this matrix and eigen_mat are the same so no copy is required)
ctx	The eigensolver context which corresponds to the MOOSE EigenProblem
tag	The MOOSE TagID for the matrix

Definition at line 679 of file SlepcSupport.C.

{
  EigenProblem * eigen_problem = static_cast<EigenProblem *>(ctx);
  auto & nl = eigen_problem->getCurrentNonlinearEigenSystem();
  auto & sys = nl.sys();
  auto & dof_map = sys.get_dof_map();
 
#ifndef NDEBUG
  auto & petsc_all_dofs_mat = cast_ref<PetscMatrix<Number> &>(all_dofs_mat);
  mooseAssert(
      !dof_map.n_constrained_dofs() == (eigen_mat == petsc_all_dofs_mat.mat()),
      "If we do not have constrained dofs, then eigen_mat and all_dofs_mat should be the same. "
      "Conversely, if we do have constrained dofs, they must be different");
#endif
 
  updateCurrentLocalSolution(sys, x);
 
  if (!eigen_problem->constJacobian())
    all_dofs_mat.zero();
 
  eigen_problem->computeJacobianTag(*sys.current_local_solution.get(), all_dofs_mat, tag);
 
  if (dof_map.n_constrained_dofs())
  {
    PetscMatrix<Number> wrapped_eigen_mat(eigen_mat, sys.comm());
    sys.copy_super_to_sub(all_dofs_mat, wrapped_eigen_mat);
  }
}

Referenced by mooseSlepcEigenFormJacobianA(), and mooseSlepcEigenFormJacobianB().

mooseSlepcEigenFormFunctionA()

PetscErrorCode Moose::SlepcSupport::mooseSlepcEigenFormFunctionA	(	SNES	snes,
		Vec	x,
		Vec	r,
		void *	ctx
	)

Form function residual Ax.

It is a SLEPc callback

Definition at line 909 of file SlepcSupport.C.

{
  PetscFunctionBegin;
 
  EigenProblem * eigen_problem = static_cast<EigenProblem *>(ctx);
  NonlinearEigenSystem & eigen_nl = eigen_problem->getCurrentNonlinearEigenSystem();
 
  if (eigen_problem->solverParams(eigen_nl.number())._eigen_matrix_vector_mult &&
      (eigen_problem->onLinearSolver() || eigen_problem->constantMatrices()))
  {
    EPS eps = eigen_nl.getEPS();
    Mat A;
    ST st;
 
    LibmeshPetscCallQ(mooseEPSFormMatrices(*eigen_problem, eps, x, ctx));
 
    // Rest ST state so that we can restrieve matrices
    LibmeshPetscCallQ(EPSGetST(eps, &st));
    LibmeshPetscCallQ(STResetMatrixState(st));
    LibmeshPetscCallQ(EPSGetOperators(eps, &A, NULL));
 
    LibmeshPetscCallQ(MatMult(A, x, r));
 
    PetscFunctionReturn(PETSC_SUCCESS);
  }
 
  moosePetscSNESFormFunction(snes, x, r, ctx, eigen_nl.nonEigenVectorTag());
 
  PetscFunctionReturn(PETSC_SUCCESS);
}

Referenced by attachCallbacksToMat().

mooseSlepcEigenFormFunctionAB()

PetscErrorCode Moose::SlepcSupport::mooseSlepcEigenFormFunctionAB	(	SNES	snes,
		Vec	x,
		Vec	Ax,
		Vec	Bx,
		void *	ctx
	)

Form function residual Ax-Bx.

It is a SLEPc callback

Definition at line 983 of file SlepcSupport.C.

{
  PetscFunctionBegin;
 
  EigenProblem * eigen_problem = static_cast<EigenProblem *>(ctx);
  NonlinearEigenSystem & eigen_nl = eigen_problem->getCurrentNonlinearEigenSystem();
  auto & sys = eigen_nl.sys();
  auto & dof_map = sys.get_dof_map();
 
  if (eigen_problem->solverParams(eigen_nl.number())._eigen_matrix_vector_mult &&
      (eigen_problem->onLinearSolver() || eigen_problem->constantMatrices()))
  {
    EPS eps = eigen_nl.getEPS();
    ST st;
    Mat A, B;
 
    LibmeshPetscCallQ(mooseEPSFormMatrices(*eigen_problem, eps, x, ctx));
 
    // Rest ST state so that we can restrieve matrices
    LibmeshPetscCallQ(EPSGetST(eps, &st));
    LibmeshPetscCallQ(STResetMatrixState(st));
 
    LibmeshPetscCallQ(EPSGetOperators(eps, &A, &B));
 
    LibmeshPetscCallQ(MatMult(A, x, Ax));
    LibmeshPetscCallQ(MatMult(B, x, Bx));
 
    if (eigen_problem->negativeSignEigenKernel())
      LibmeshPetscCallQ(VecScale(Bx, -1.));
 
    if (eigen_problem->bxNormProvided())
    {
      // User has provided a postprocessor. We need it updated
      updateCurrentLocalSolution(sys, x);
      eigen_problem->execute(EXEC_LINEAR);
    }
 
    PetscFunctionReturn(PETSC_SUCCESS);
  }
 
  updateCurrentLocalSolution(sys, x);
  auto AX = createWrappedResidual(sys, Ax);
  auto BX = createWrappedResidual(sys, Bx);
 
  eigen_problem->computeResidualAB(*sys.current_local_solution.get(),
                                   *AX,
                                   *BX,
                                   eigen_nl.nonEigenVectorTag(),
                                   eigen_nl.eigenVectorTag());
 
  AX->close();
  BX->close();
 
  if (dof_map.n_constrained_dofs())
  {
    PetscVector<Number> sub_Ax(Ax, sys.comm());
    sys.copy_super_to_sub(*AX, sub_Ax);
    PetscVector<Number> sub_Bx(Bx, sys.comm());
    sys.copy_super_to_sub(*BX, sub_Bx);
  }
 
  if (eigen_problem->negativeSignEigenKernel())
    LibmeshPetscCallQ(VecScale(Bx, -1.));
 
  PetscFunctionReturn(PETSC_SUCCESS);
}

Referenced by attachCallbacksToMat().

mooseSlepcEigenFormFunctionB()

PetscErrorCode Moose::SlepcSupport::mooseSlepcEigenFormFunctionB	(	SNES	snes,
		Vec	x,
		Vec	r,
		void *	ctx
	)

Form function residual Bx.

It is a SLEPc callback

Definition at line 941 of file SlepcSupport.C.

{
  PetscFunctionBegin;
 
  EigenProblem * eigen_problem = static_cast<EigenProblem *>(ctx);
  NonlinearEigenSystem & eigen_nl = eigen_problem->getCurrentNonlinearEigenSystem();
 
  if (eigen_problem->solverParams(eigen_nl.number())._eigen_matrix_vector_mult &&
      (eigen_problem->onLinearSolver() || eigen_problem->constantMatrices()))
  {
    EPS eps = eigen_nl.getEPS();
    ST st;
    Mat B;
 
    LibmeshPetscCallQ(mooseEPSFormMatrices(*eigen_problem, eps, x, ctx));
 
    // Rest ST state so that we can restrieve matrices
    LibmeshPetscCallQ(EPSGetST(eps, &st));
    LibmeshPetscCallQ(STResetMatrixState(st));
    LibmeshPetscCallQ(EPSGetOperators(eps, NULL, &B));
 
    LibmeshPetscCallQ(MatMult(B, x, r));
 
    if (eigen_problem->bxNormProvided())
    {
      // User has provided a postprocessor. We need it updated
      updateCurrentLocalSolution(eigen_nl.sys(), x);
      eigen_problem->execute(EXEC_LINEAR);
    }
  }
  else
    moosePetscSNESFormFunction(snes, x, r, ctx, eigen_nl.eigenVectorTag());
 
  if (eigen_problem->negativeSignEigenKernel())
  {
    LibmeshPetscCallQ(VecScale(r, -1.));
  }
 
  PetscFunctionReturn(PETSC_SUCCESS);
}

Referenced by attachCallbacksToMat().

mooseSlepcEigenFormFunctionMFFD()

PetscErrorCode Moose::SlepcSupport::mooseSlepcEigenFormFunctionMFFD	(	void *	ctx,
		Vec	x,
		Vec	r
	)

Function call for MFFD.

Definition at line 747 of file SlepcSupport.C.

{
  PetscErrorCode (*func)(SNES, Vec, Vec, void *);
  void * fctx;
  PetscFunctionBegin;
 
  EigenProblem * eigen_problem = static_cast<EigenProblem *>(ctx);
  NonlinearEigenSystem & eigen_nl = eigen_problem->getCurrentNonlinearEigenSystem();
  SNES snes = eigen_nl.getSNES();
 
  eigen_problem->onLinearSolver(true);
 
  LibmeshPetscCallQ(SNESGetFunction(snes, NULL, &func, &fctx));
  if (fctx != ctx)
  {
    SETERRQ(
        PetscObjectComm((PetscObject)snes), PETSC_ERR_ARG_INCOMP, "Contexts are not consistent \n");
  }
  LibmeshPetscCallQ((*func)(snes, x, r, ctx));
 
  eigen_problem->onLinearSolver(false);
 
  PetscFunctionReturn(PETSC_SUCCESS);
}

Referenced by mooseSlepcEigenFormJacobianA().

mooseSlepcEigenFormJacobianA()

PetscErrorCode Moose::SlepcSupport::mooseSlepcEigenFormJacobianA	(	SNES	snes,
		Vec	x,
		Mat	jac,
		Mat	pc,
		void *	ctx
	)

Form Jacobian matrix A.

It is a SLEPc callback

Definition at line 773 of file SlepcSupport.C.

{
  PetscBool jisshell, pisshell;
  PetscBool jismffd;
 
  PetscFunctionBegin;
 
  EigenProblem * eigen_problem = static_cast<EigenProblem *>(ctx);
  NonlinearEigenSystem & eigen_nl = eigen_problem->getCurrentNonlinearEigenSystem();
  auto & sys = eigen_nl.sys();
 
  // If both jacobian and preconditioning are shell matrices,
  // and then assemble them and return
  LibmeshPetscCallQ(PetscObjectTypeCompare((PetscObject)jac, MATSHELL, &jisshell));
  LibmeshPetscCallQ(PetscObjectTypeCompare((PetscObject)jac, MATMFFD, &jismffd));
 
  if (jismffd && eigen_problem->solverParams(eigen_nl.number())._eigen_matrix_vector_mult)
  {
    LibmeshPetscCallQ(
        MatMFFDSetFunction(jac, Moose::SlepcSupport::mooseSlepcEigenFormFunctionMFFD, ctx));
 
    EPS eps = eigen_nl.getEPS();
 
    LibmeshPetscCallQ(mooseEPSFormMatrices(*eigen_problem, eps, x, ctx));
 
    if (pc != jac)
    {
      LibmeshPetscCallQ(MatAssemblyBegin(jac, MAT_FINAL_ASSEMBLY));
      LibmeshPetscCallQ(MatAssemblyEnd(jac, MAT_FINAL_ASSEMBLY));
    }
    PetscFunctionReturn(PETSC_SUCCESS);
  }
 
  LibmeshPetscCallQ(PetscObjectTypeCompare((PetscObject)pc, MATSHELL, &pisshell));
  if ((jisshell || jismffd) && pisshell)
  {
    // Just assemble matrices and return
    LibmeshPetscCallQ(MatAssemblyBegin(jac, MAT_FINAL_ASSEMBLY));
    LibmeshPetscCallQ(MatAssemblyBegin(pc, MAT_FINAL_ASSEMBLY));
    LibmeshPetscCallQ(MatAssemblyEnd(jac, MAT_FINAL_ASSEMBLY));
    LibmeshPetscCallQ(MatAssemblyEnd(pc, MAT_FINAL_ASSEMBLY));
 
    PetscFunctionReturn(PETSC_SUCCESS);
  }
 
  // Jacobian and precond matrix are the same
  if (jac == pc)
  {
    if (!pisshell)
      moosePetscSNESFormMatrixTag(
          snes, x, pc, sys.get_matrix_A(), ctx, eigen_nl.precondMatrixTag());
 
    PetscFunctionReturn(PETSC_SUCCESS);
  }
  else
  {
    if (!jisshell && !jismffd && !pisshell) // We need to form both Jacobian and precond matrix
    {
      std::vector<Mat> mats = {jac, pc};
      std::vector<SparseMatrix<Number> *> libmesh_mats = {&sys.get_matrix_A(),
                                                          &sys.get_precond_matrix()};
      std::set<TagID> tags = {eigen_nl.nonEigenMatrixTag(), eigen_nl.precondMatrixTag()};
      moosePetscSNESFormMatricesTags(snes, x, mats, libmesh_mats, ctx, tags);
      PetscFunctionReturn(PETSC_SUCCESS);
    }
    if (!pisshell) // We need to form only precond matrix
    {
      moosePetscSNESFormMatrixTag(
          snes, x, pc, sys.get_precond_matrix(), ctx, eigen_nl.precondMatrixTag());
      LibmeshPetscCallQ(MatAssemblyBegin(jac, MAT_FINAL_ASSEMBLY));
      LibmeshPetscCallQ(MatAssemblyEnd(jac, MAT_FINAL_ASSEMBLY));
      PetscFunctionReturn(PETSC_SUCCESS);
    }
    if (!jisshell && !jismffd) // We need to form only Jacobian matrix
    {
      moosePetscSNESFormMatrixTag(
          snes, x, jac, sys.get_matrix_A(), ctx, eigen_nl.nonEigenMatrixTag());
      LibmeshPetscCallQ(MatAssemblyBegin(pc, MAT_FINAL_ASSEMBLY));
      LibmeshPetscCallQ(MatAssemblyEnd(pc, MAT_FINAL_ASSEMBLY));
      PetscFunctionReturn(PETSC_SUCCESS);
    }
  }
  PetscFunctionReturn(PETSC_SUCCESS);
}

Referenced by attachCallbacksToMat().

mooseSlepcEigenFormJacobianB()

PetscErrorCode Moose::SlepcSupport::mooseSlepcEigenFormJacobianB	(	SNES	snes,
		Vec	x,
		Mat	jac,
		Mat	pc,
		void *	ctx
	)

Form Jacobian matrix B.

It is a SLEPc callback

Definition at line 859 of file SlepcSupport.C.

{
  PetscBool jshell, pshell;
  PetscBool jismffd;
 
  PetscFunctionBegin;
 
  EigenProblem * eigen_problem = static_cast<EigenProblem *>(ctx);
  NonlinearEigenSystem & eigen_nl = eigen_problem->getCurrentNonlinearEigenSystem();
  auto & sys = eigen_nl.sys();
 
  // If both jacobian and preconditioning are shell matrices,
  // and then assemble them and return
  LibmeshPetscCallQ(PetscObjectTypeCompare((PetscObject)jac, MATSHELL, &jshell));
  LibmeshPetscCallQ(PetscObjectTypeCompare((PetscObject)jac, MATMFFD, &jismffd));
  LibmeshPetscCallQ(PetscObjectTypeCompare((PetscObject)pc, MATSHELL, &pshell));
  if ((jshell || jismffd) && pshell)
  {
    // Just assemble matrices and return
    LibmeshPetscCallQ(MatAssemblyBegin(jac, MAT_FINAL_ASSEMBLY));
    LibmeshPetscCallQ(MatAssemblyBegin(pc, MAT_FINAL_ASSEMBLY));
    LibmeshPetscCallQ(MatAssemblyEnd(jac, MAT_FINAL_ASSEMBLY));
    LibmeshPetscCallQ(MatAssemblyEnd(pc, MAT_FINAL_ASSEMBLY));
 
    PetscFunctionReturn(PETSC_SUCCESS);
  }
 
  if (jac != pc && (!jshell && !jshell))
    SETERRQ(PetscObjectComm((PetscObject)snes),
            PETSC_ERR_ARG_INCOMP,
            "Jacobian and precond matrices should be the same for eigen kernels \n");
 
  moosePetscSNESFormMatrixTag(snes, x, pc, sys.get_matrix_B(), ctx, eigen_nl.eigenMatrixTag());
 
  if (eigen_problem->negativeSignEigenKernel())
  {
    LibmeshPetscCallQ(MatScale(pc, -1.));
  }
 
  PetscFunctionReturn(PETSC_SUCCESS);
}

Referenced by attachCallbacksToMat().

mooseSlepcEigenFormNorm()

PetscErrorCode Moose::SlepcSupport::mooseSlepcEigenFormNorm	(	SNES	,
		Vec	,
		PetscReal *	norm,
		void *	ctx
	)

Definition at line 1051 of file SlepcSupport.C.

{
  PetscFunctionBegin;
  auto * const eigen_problem = static_cast<EigenProblem *>(ctx);
  *norm = eigen_problem->formNorm();
  PetscFunctionReturn(PETSC_SUCCESS);
}

Referenced by attachCallbacksToMat().

mooseSlepcEPSGetSNES()

PetscErrorCode Moose::SlepcSupport::mooseSlepcEPSGetSNES	(	EPS	eps,
		SNES *	snes
	)

Retrieve SNES from EPS.

It makes sense only for Newton eigenvalue solvers

Definition at line 1309 of file SlepcSupport.C.

{
  PetscBool same, nonlinear;
 
  PetscFunctionBegin;
  LibmeshPetscCallQ(PetscObjectTypeCompare((PetscObject)eps, EPSPOWER, &same));
 
  if (!same)
    mooseError("It is not eps power, and there is no snes");
 
  LibmeshPetscCallQ(EPSPowerGetNonlinear(eps, &nonlinear));
 
  if (!nonlinear)
    mooseError("It is not a nonlinear eigen solver");
 
  LibmeshPetscCallQ(EPSPowerGetSNES(eps, snes));
 
  PetscFunctionReturn(PETSC_SUCCESS);
}

Referenced by NonlinearEigenSystem::getSNES(), mooseSlepcEPSSNESKSPSetPCSide(), mooseSlepcEPSSNESSetCustomizePC(), and mooseSlepcEPSSNESSetUpOptionPrefix().

mooseSlepcEPSMonitor()

PetscErrorCode Moose::SlepcSupport::mooseSlepcEPSMonitor	(	EPS	eps,
		PetscInt	its,
		PetscInt	nconv,
		PetscScalar *	eigr,
		PetscScalar *	eigi,
		PetscReal *	errest,
		PetscInt	nest,
		void *	mctx
	)

A customized solver monitor to print out eigenvalue.

Definition at line 1385 of file SlepcSupport.C.

{
  ST st;
  PetscScalar eigenr, eigeni;
 
  PetscFunctionBegin;
  EigenProblem * eigen_problem = static_cast<EigenProblem *>(mctx);
  auto & console = eigen_problem->console();
 
  auto inverse = eigen_problem->outputInverseEigenvalue();
  LibmeshPetscCallQ(EPSGetST(eps, &st));
  eigenr = eigr[0];
  eigeni = eigi[0];
  // Make the eigenvalue consistent with shift type
  LibmeshPetscCallQ(STBackTransform(st, 1, &eigenr, &eigeni));
 
  auto eigenvalue = inverse ? 1.0 / eigenr : eigenr;
 
  // The term "k-eigenvalue" is adopted from the neutronics community.
  console << " Iteration " << its << std::setprecision(10) << std::fixed
          << (inverse ? " k-eigenvalue = " : " eigenvalue = ") << eigenvalue << std::endl;
 
  PetscFunctionReturn(PETSC_SUCCESS);
}

Referenced by SlepcEigenSolverConfiguration::configure_solver().

mooseSlepcEPSSNESKSPSetPCSide()

PetscErrorCode Moose::SlepcSupport::mooseSlepcEPSSNESKSPSetPCSide	(	FEProblemBase &	problem,
		EPS	eps
	)

Allow users to specify PC side.

By default, we apply PC from the right side. It is consistent with the nonlinear solver.

Definition at line 1367 of file SlepcSupport.C.

{
  SNES snes;
  KSP ksp;
 
  PetscFunctionBegin;
  // Get SNES from EPS
  LibmeshPetscCallQ(mooseSlepcEPSGetSNES(eps, &snes));
  // Get KSP from SNES
  LibmeshPetscCallQ(SNESGetKSP(snes, &ksp));
 
  Moose::PetscSupport::petscSetDefaultPCSide(problem, ksp);
 
  Moose::PetscSupport::petscSetDefaultKSPNormType(problem, ksp);
  PetscFunctionReturn(PETSC_SUCCESS);
}

Referenced by SlepcEigenSolverConfiguration::configure_solver().

mooseSlepcEPSSNESSetCustomizePC()

PetscErrorCode Moose::SlepcSupport::mooseSlepcEPSSNESSetCustomizePC

(

EPS

eps

)

Attach a customized PC.

It is useful when you want to use PBP or other customized preconditioners

Definition at line 1348 of file SlepcSupport.C.

{
  SNES snes;
  KSP ksp;
  PC pc;
 
  PetscFunctionBegin;
  // Get SNES from EPS
  LibmeshPetscCallQ(mooseSlepcEPSGetSNES(eps, &snes));
  // Get KSP from SNES
  LibmeshPetscCallQ(SNESGetKSP(snes, &ksp));
  // Get PC from KSP
  LibmeshPetscCallQ(KSPGetPC(ksp, &pc));
  // Set PC type
  LibmeshPetscCallQ(PCSetType(pc, "moosepc"));
  PetscFunctionReturn(PETSC_SUCCESS);
}

Referenced by SlepcEigenSolverConfiguration::configure_solver().

mooseSlepcEPSSNESSetUpOptionPrefix()

PetscErrorCode Moose::SlepcSupport::mooseSlepcEPSSNESSetUpOptionPrefix

(

EPS

eps

)

Get rid of prefix "-eps_power" for SNES, KSP, PC, etc.

Definition at line 1330 of file SlepcSupport.C.

{
  SNES snes;
  const char * prefix = nullptr;
 
  PetscFunctionBegin;
  LibmeshPetscCallQ(mooseSlepcEPSGetSNES(eps, &snes));
  // There is an extra "eps_power" in snes that users do not like it.
  // Let us remove that from snes.
  // Retrieve option prefix from EPS
  LibmeshPetscCallQ(PetscObjectGetOptionsPrefix((PetscObject)eps, &prefix));
  // Set option prefix to SNES
  LibmeshPetscCallQ(SNESSetOptionsPrefix(snes, prefix));
 
  PetscFunctionReturn(PETSC_SUCCESS);
}

Referenced by SlepcEigenSolverConfiguration::configure_solver().

mooseSlepcStoppingTest()

PetscErrorCode Moose::SlepcSupport::mooseSlepcStoppingTest	(	EPS	eps,
		PetscInt	its,
		PetscInt	max_it,
		PetscInt	nconv,
		PetscInt	nev,
		EPSConvergedReason *	reason,
		void *	ctx
	)

A customized convergence checker.

We need to make solver as converged when doing free power iteration.

Definition at line 1281 of file SlepcSupport.C.

{
  EigenProblem * eigen_problem = static_cast<EigenProblem *>(ctx);
 
  PetscFunctionBegin;
  LibmeshPetscCallQ(EPSStoppingBasic(eps, its, max_it, nconv, nev, reason, NULL));
 
  // If we do free power iteration, we need to mark the solver as converged.
  // It is because SLEPc does not offer a way to copy unconverged solution.
  // If the solver is not marked as "converged", we have no way to get solution
  // from slepc. Note marking as "converged" has no side-effects at all for us.
  // If free power iteration is used as a stand-alone solver, we won't trigger
  // as "doFreePowerIteration()" is false.
  if (eigen_problem->doFreePowerIteration() && its == max_it && *reason <= 0)
  {
    *reason = EPS_CONVERGED_USER;
    eps->nconv = 1;
  }
  PetscFunctionReturn(PETSC_SUCCESS);
}

Referenced by SlepcEigenSolverConfiguration::configure_solver().

PCApply_MoosePC()

PetscErrorCode Moose::SlepcSupport::PCApply_MoosePC	(	PC	pc,
		Vec	x,
		Vec	y
	)

Preconditioner application.

It call libMesh preconditioner to implement this.

Definition at line 1219 of file SlepcSupport.C.

{
  void * ctx;
  Mat Amat, Pmat;
  PetscContainer container;
 
  PetscFunctionBegin;
  LibmeshPetscCallQ(PCGetOperators(pc, &Amat, &Pmat));
  LibmeshPetscCallQ(
      PetscObjectQuery((PetscObject)Pmat, "formFunctionCtx", (PetscObject *)&container));
  if (container)
    LibmeshPetscCallQ(PetscContainerGetPointer(container, &ctx));
  else
    mooseError(" Can not find a context \n");
 
  EigenProblem * eigen_problem = static_cast<EigenProblem *>(ctx);
  NonlinearEigenSystem & nl_eigen = eigen_problem->getCurrentNonlinearEigenSystem();
  auto preconditioner = nl_eigen.preconditioner();
 
  if (!preconditioner)
    mooseError("There is no moose preconditioner in nonlinear eigen system \n");
 
  PetscVector<Number> x_vec(x, preconditioner->comm());
  PetscVector<Number> y_vec(y, preconditioner->comm());
 
  preconditioner->apply(x_vec, y_vec);
 
  PetscFunctionReturn(PETSC_SUCCESS);
}

Referenced by PCCreate_MoosePC().

PCCreate_MoosePC()

PETSC_EXTERN PetscErrorCode Moose::SlepcSupport::PCCreate_MoosePC

(

PC

pc

)

Create a preconditioner from moose side.

It is used to attach moose preconditioner

Definition at line 1184 of file SlepcSupport.C.

{
  PetscFunctionBegin;
 
  pc->ops->view = PCView_MoosePC;
  pc->ops->destroy = PCDestroy_MoosePC;
  pc->ops->setup = PCSetUp_MoosePC;
  pc->ops->apply = PCApply_MoosePC;
 
  PetscFunctionReturn(PETSC_SUCCESS);
}

Referenced by registerPCToPETSc().

PCDestroy_MoosePC()

PetscErrorCode Moose::SlepcSupport::PCDestroy_MoosePC

(

PC

pc

)

Destroy preconditioner.

Definition at line 1197 of file SlepcSupport.C.

{
  PetscFunctionBegin;
  /* We do not need to do anything right now, but later we may have some data we need to free here
   */
  PetscFunctionReturn(PETSC_SUCCESS);
}

Referenced by PCCreate_MoosePC().

PCSetUp_MoosePC()

PetscErrorCode Moose::SlepcSupport::PCSetUp_MoosePC

(

PC

pc

)

Setup preconditioner.

Definition at line 1250 of file SlepcSupport.C.

{
  void * ctx;
  Mat Amat, Pmat;
  PetscContainer container;
 
  PetscFunctionBegin;
  LibmeshPetscCallQ(PCGetOperators(pc, &Amat, &Pmat));
  LibmeshPetscCallQ(
      PetscObjectQuery((PetscObject)Pmat, "formFunctionCtx", (PetscObject *)&container));
  if (container)
    LibmeshPetscCallQ(PetscContainerGetPointer(container, &ctx));
  else
    mooseError(" Can not find a context \n");
 
  EigenProblem * eigen_problem = static_cast<EigenProblem *>(ctx);
  NonlinearEigenSystem & nl_eigen = eigen_problem->getCurrentNonlinearEigenSystem();
  Preconditioner<Number> * preconditioner = nl_eigen.preconditioner();
 
  if (!preconditioner)
    mooseError("There is no moose preconditioner in nonlinear eigen system \n");
 
  if (!preconditioner->initialized())
    preconditioner->init();
 
  preconditioner->setup();
 
  PetscFunctionReturn(PETSC_SUCCESS);
}

Referenced by PCCreate_MoosePC().

PCView_MoosePC()

PetscErrorCode Moose::SlepcSupport::PCView_MoosePC	(	PC	pc,
		PetscViewer	viewer
	)

View preconditioner.

Definition at line 1206 of file SlepcSupport.C.

{
  PetscBool iascii;
 
  PetscFunctionBegin;
  LibmeshPetscCallQ(PetscObjectTypeCompare((PetscObject)viewer, PETSCVIEWERASCII, &iascii));
  if (iascii)
    LibmeshPetscCallQ(PetscViewerASCIIPrintf(viewer, "  %s\n", "moosepc"));
 
  PetscFunctionReturn(PETSC_SUCCESS);
}

Referenced by PCCreate_MoosePC().

registerPCToPETSc()

PETSC_EXTERN PetscErrorCode Moose::SlepcSupport::registerPCToPETSc

(

)

Let PETSc know there is a preconditioner.

Definition at line 1174 of file SlepcSupport.C.

{
  PetscFunctionBegin;
 
  LibmeshPetscCallQ(PCRegister("moosepc", PCCreate_MoosePC));
 
  PetscFunctionReturn(PETSC_SUCCESS);
}

Referenced by NonlinearEigenSystem::attachPreconditioner().

setEigenProblemOptions()

void Moose::SlepcSupport::setEigenProblemOptions	(	SolverParams &	solver_params,
		const MultiMooseEnum &	dont_add_these_options
	)

Definition at line 296 of file SlepcSupport.C.

{
  switch (solver_params._eigen_problem_type)
  {
    case Moose::EPT_HERMITIAN:
      Moose::PetscSupport::setSinglePetscOptionIfAppropriate(dont_add_these_options,
                                                             "-eps_hermitian");
      break;
 
    case Moose::EPT_NON_HERMITIAN:
      Moose::PetscSupport::setSinglePetscOptionIfAppropriate(dont_add_these_options,
                                                             "-eps_non_hermitian");
      break;
 
    case Moose::EPT_GEN_HERMITIAN:
      Moose::PetscSupport::setSinglePetscOptionIfAppropriate(dont_add_these_options,
                                                             "-eps_gen_hermitian");
      break;
 
    case Moose::EPT_GEN_INDEFINITE:
      Moose::PetscSupport::setSinglePetscOptionIfAppropriate(dont_add_these_options,
                                                             "-eps_gen_indefinite");
      break;
 
    case Moose::EPT_GEN_NON_HERMITIAN:
      Moose::PetscSupport::setSinglePetscOptionIfAppropriate(dont_add_these_options,
                                                             "-eps_gen_non_hermitian");
      break;
 
    case Moose::EPT_POS_GEN_NON_HERMITIAN:
      Moose::PetscSupport::setSinglePetscOptionIfAppropriate(dont_add_these_options,
                                                             "-eps_pos_gen_non_hermitian");
      break;
 
    case Moose::EPT_SLEPC_DEFAULT:
      break;
 
    default:
      mooseError("Unknown eigen solver type \n");
  }
}

Referenced by slepcSetOptions().

setEigenProblemSolverParams()

void Moose::SlepcSupport::setEigenProblemSolverParams	(	EigenProblem &	eigen_problem,
		const InputParameters &	params
	)

Retrieve eigen problem params from 'params', and then set these params into SolverParams.

Definition at line 210 of file SlepcSupport.C.

{
  for (const auto i : make_range(eigen_problem.numNonlinearSystems()))
  {
    const std::string & eigen_problem_type = params.get<MooseEnum>("eigen_problem_type");
    if (!eigen_problem_type.empty())
      eigen_problem.solverParams(i)._eigen_problem_type =
          Moose::stringToEnum<Moose::EigenProblemType>(eigen_problem_type);
    else
      mooseError("Have to specify a valid eigen problem type");
 
    const std::string & which_eigen_pairs = params.get<MooseEnum>("which_eigen_pairs");
    if (!which_eigen_pairs.empty())
      eigen_problem.solverParams(i)._which_eigen_pairs =
          Moose::stringToEnum<Moose::WhichEigenPairs>(which_eigen_pairs);
 
    // Set necessary parameters used in EigenSystem::solve(),
    // i.e. the number of requested eigenpairs nev and the number
    // of basis vectors ncv used in the solution algorithm. Note that
    // ncv >= nev must hold and ncv >= 2*nev is recommended
    unsigned int n_eigen_pairs = params.get<unsigned int>("n_eigen_pairs");
    unsigned int n_basis_vectors = params.get<unsigned int>("n_basis_vectors");
 
    eigen_problem.setNEigenPairsRequired(n_eigen_pairs);
 
    eigen_problem.es().parameters.set<unsigned int>("eigenpairs") = n_eigen_pairs;
 
    // If the subspace dimension is too small, we increase it automatically
    if (subspace_factor * n_eigen_pairs > n_basis_vectors)
    {
      n_basis_vectors = subspace_factor * n_eigen_pairs;
      mooseWarning(
          "Number of subspaces in Eigensolver is changed by moose because the value you set "
          "is too small");
    }
 
    eigen_problem.es().parameters.set<unsigned int>("basis vectors") = n_basis_vectors;
 
    // Operators A and B are formed as shell matrices
    eigen_problem.solverParams(i)._eigen_matrix_free = params.get<bool>("matrix_free");
 
    // Preconditioning is formed as a shell matrix
    eigen_problem.solverParams(i)._precond_matrix_free = params.get<bool>("precond_matrix_free");
 
    if (params.get<MooseEnum>("solve_type") == "PJFNK")
    {
      eigen_problem.solverParams(i)._eigen_matrix_free = true;
    }
    if (params.get<MooseEnum>("solve_type") == "JFNK")
    {
      eigen_problem.solverParams(i)._eigen_matrix_free = true;
      eigen_problem.solverParams(i)._precond_matrix_free = true;
    }
    // We need matrices so that we can implement residual evaluations
    if (params.get<MooseEnum>("solve_type") == "PJFNKMO")
    {
      eigen_problem.solverParams(i)._eigen_matrix_free = true;
      eigen_problem.solverParams(i)._precond_matrix_free = false;
      eigen_problem.solverParams(i)._eigen_matrix_vector_mult = true;
      // By default, we need to form full matrices, otherwise residual
      // evaluations will not be accurate
      eigen_problem.setCoupling(Moose::COUPLING_FULL);
    }
  }
 
  eigen_problem.constantMatrices(params.get<bool>("constant_matrices"));
 
  if (eigen_problem.constantMatrices() && params.get<MooseEnum>("solve_type") != "PJFNKMO")
  {
    mooseError("constant_matrices flag is only valid for solve type: PJFNKMO");
  }
}

Referenced by EigenProblemSolve::EigenProblemSolve().

setEigenSolverOptions()

void Moose::SlepcSupport::setEigenSolverOptions	(	SolverParams &	solver_params,
		const InputParameters &	params
	)

Definition at line 492 of file SlepcSupport.C.

{
  // Avoid unused variable warnings when you have SLEPc but not PETSc-dev.
  libmesh_ignore(params);
 
  switch (solver_params._eigen_solve_type)
  {
    case Moose::EST_POWER:
      Moose::PetscSupport::setSinglePetscOption("-eps_type", "power");
      break;
 
    case Moose::EST_ARNOLDI:
      Moose::PetscSupport::setSinglePetscOption("-eps_type", "arnoldi");
      break;
 
    case Moose::EST_KRYLOVSCHUR:
      Moose::PetscSupport::setSinglePetscOption("-eps_type", "krylovschur");
      break;
 
    case Moose::EST_JACOBI_DAVIDSON:
      Moose::PetscSupport::setSinglePetscOption("-eps_type", "jd");
      break;
 
    case Moose::EST_NONLINEAR_POWER:
      setNonlinearPowerOptions(solver_params);
      break;
 
    case Moose::EST_NEWTON:
      setNewtonPetscOptions(solver_params, params);
      break;
 
    case Moose::EST_PJFNK:
      solver_params._eigen_matrix_free = true;
      solver_params._customized_pc_for_eigen = false;
      setNewtonPetscOptions(solver_params, params);
      break;
 
    case Moose::EST_JFNK:
      solver_params._eigen_matrix_free = true;
      solver_params._customized_pc_for_eigen = true;
      setNewtonPetscOptions(solver_params, params);
      break;
 
    case Moose::EST_PJFNKMO:
      solver_params._eigen_matrix_free = true;
      solver_params._customized_pc_for_eigen = false;
      solver_params._eigen_matrix_vector_mult = true;
      setNewtonPetscOptions(solver_params, params);
      break;
 
    default:
      mooseError("Unknown eigen solver type \n");
  }
}

Referenced by slepcSetOptions().

setFreeNonlinearPowerIterations()

void Moose::SlepcSupport::setFreeNonlinearPowerIterations

(

unsigned int

free_power_iterations

)

Set SLEPc/PETSc options to trigger free power iteration.

Definition at line 402 of file SlepcSupport.C.

{
  Moose::PetscSupport::setSinglePetscOption("-eps_power_update", "0");
  Moose::PetscSupport::setSinglePetscOption("-snes_max_it", "2");
  // During each power iteration, we want solver converged unless linear solver does not
  // work. We here use a really loose tolerance for this purpose.
  // -snes_no_convergence_test is a perfect option, but it was removed from PETSc
  Moose::PetscSupport::setSinglePetscOption("-snes_rtol", "0.99999999999");
  Moose::PetscSupport::setSinglePetscOption("-eps_max_it", stringify(free_power_iterations));
  // We always want the number of free power iterations respected so we don't want to stop early if
  // we've satisfied a convergence criterion. Consequently we make this tolerance very tight
  Moose::PetscSupport::setSinglePetscOption("-eps_tol", "1e-50");
}

Referenced by EigenProblem::doFreeNonlinearPowerIterations().

setNewtonPetscOptions()

void Moose::SlepcSupport::setNewtonPetscOptions	(	SolverParams &	solver_params,
		const InputParameters &	params
	)

Definition at line 429 of file SlepcSupport.C.

{
#if !SLEPC_VERSION_LESS_THAN(3, 8, 0) || !PETSC_VERSION_RELEASE
  // Whether or not we need to involve an initial inverse power
  bool initial_power = params.get<bool>("_newton_inverse_power");
 
  Moose::PetscSupport::setSinglePetscOption("-eps_type", "power");
  Moose::PetscSupport::setSinglePetscOption("-eps_power_nonlinear", "1");
  Moose::PetscSupport::setSinglePetscOption("-eps_power_update", "1");
  // Only one outer iteration in EPS is allowed when Newton/PJFNK/JFNK
  // is used as the eigen solver
  Moose::PetscSupport::setSinglePetscOption("-eps_max_it", "1");
  if (initial_power)
  {
    Moose::PetscSupport::setSinglePetscOption("-init_eps_power_snes_max_it", "1");
    Moose::PetscSupport::setSinglePetscOption("-init_eps_power_ksp_rtol", "1e-2");
    Moose::PetscSupport::setSinglePetscOption(
        "-init_eps_max_it", stringify(params.get<unsigned int>("free_power_iterations")));
  }
  Moose::PetscSupport::setSinglePetscOption("-eps_target_magnitude", "");
  if (solver_params._eigen_matrix_free)
  {
    Moose::PetscSupport::setSinglePetscOption("-snes_mf_operator", "1");
    if (initial_power)
      Moose::PetscSupport::setSinglePetscOption("-init_eps_power_snes_mf_operator", "1");
  }
  else
  {
    Moose::PetscSupport::setSinglePetscOption("-snes_mf_operator", "0");
    if (initial_power)
      Moose::PetscSupport::setSinglePetscOption("-init_eps_power_snes_mf_operator", "0");
  }
#if PETSC_RELEASE_LESS_THAN(3, 13, 0)
  Moose::PetscSupport::setSinglePetscOption("-st_type", "sinvert");
  if (initial_power)
    Moose::PetscSupport::setSinglePetscOption("-init_st_type", "sinvert");
#endif
#else
  mooseError("Newton-based eigenvalue solver requires SLEPc 3.7.3 or higher");
#endif
}

Referenced by setEigenSolverOptions().

setNonlinearPowerOptions()

void Moose::SlepcSupport::setNonlinearPowerOptions

(

SolverParams &

solver_params

)

Definition at line 472 of file SlepcSupport.C.

{
#if !SLEPC_VERSION_LESS_THAN(3, 8, 0) || !PETSC_VERSION_RELEASE
  Moose::PetscSupport::setSinglePetscOption("-eps_type", "power");
  Moose::PetscSupport::setSinglePetscOption("-eps_power_nonlinear", "1");
  Moose::PetscSupport::setSinglePetscOption("-eps_target_magnitude", "");
  if (solver_params._eigen_matrix_free)
    Moose::PetscSupport::setSinglePetscOption("-snes_mf_operator", "1");
  else
    Moose::PetscSupport::setSinglePetscOption("-snes_mf_operator", "0");
 
#if PETSC_RELEASE_LESS_THAN(3, 13, 0)
  Moose::PetscSupport::setSinglePetscOption("-st_type", "sinvert");
#endif
#else
  mooseError("Nonlinear Inverse Power requires SLEPc 3.7.3 or higher");
#endif
}

Referenced by setEigenSolverOptions().

setOperationsForShellMat()

void Moose::SlepcSupport::setOperationsForShellMat	(	EigenProblem &	eigen_problem,
		Mat	mat,
		bool	eigen
	)

Set operations to shell mat.

Definition at line 1163 of file SlepcSupport.C.

{
  LibmeshPetscCallA(eigen_problem.comm().get(), MatShellSetContext(mat, &eigen_problem));
  LibmeshPetscCallA(eigen_problem.comm().get(),
                    MatShellSetOperation(mat,
                                         MATOP_MULT,
                                         eigen ? (void (*)(void))mooseMatMult_Eigen
                                               : (void (*)(void))mooseMatMult_NonEigen));
}

Referenced by NonlinearEigenSystem::attachSLEPcCallbacks().

setSlepcEigenSolverTolerances()

void Moose::SlepcSupport::setSlepcEigenSolverTolerances	(	EigenProblem &	eigen_problem,
		const SolverParams &	solver_params,
		const InputParameters &	params
	)

Control eigen solver tolerances via SLEPc options.

Definition at line 126 of file SlepcSupport.C.

{
  const auto & dont_add_these_options = eigen_problem.getPetscOptions().dont_add_these_options;
  const auto prefix_with_dash = '-' + solver_params._prefix;
 
  mooseAssert(solver_params._solver_sys_num != libMesh::invalid_uint,
              "The solver system number must be initialized");
 
  Moose::PetscSupport::setSinglePetscOptionIfAppropriate(dont_add_these_options,
                                                         prefix_with_dash + "eps_tol",
                                                         stringify(params.get<Real>("eigen_tol")));
 
  Moose::PetscSupport::setSinglePetscOptionIfAppropriate(
      dont_add_these_options,
      prefix_with_dash + "eps_max_it",
      stringify(params.get<unsigned int>("eigen_max_its")));
 
  // if it is a nonlinear eigenvalue solver, we need to set tolerances for nonlinear solver and
  // linear solver
  if (eigen_problem.isNonlinearEigenvalueSolver(solver_params._solver_sys_num))
  {
    // nonlinear solver tolerances
    Moose::PetscSupport::setSinglePetscOptionIfAppropriate(
        dont_add_these_options,
        prefix_with_dash + "snes_max_it",
        stringify(params.get<unsigned int>("nl_max_its")));
 
    Moose::PetscSupport::setSinglePetscOptionIfAppropriate(
        dont_add_these_options,
        prefix_with_dash + "snes_max_funcs",
        stringify(params.get<unsigned int>("nl_max_funcs")));
 
    Moose::PetscSupport::setSinglePetscOptionIfAppropriate(
        dont_add_these_options,
        prefix_with_dash + "snes_atol",
        stringify(params.get<Real>("nl_abs_tol")));
 
    Moose::PetscSupport::setSinglePetscOptionIfAppropriate(
        dont_add_these_options,
        prefix_with_dash + "snes_rtol",
        stringify(params.get<Real>("nl_rel_tol")));
 
    Moose::PetscSupport::setSinglePetscOptionIfAppropriate(
        dont_add_these_options,
        prefix_with_dash + "snes_stol",
        stringify(params.get<Real>("nl_rel_step_tol")));
 
    // linear solver
    Moose::PetscSupport::setSinglePetscOptionIfAppropriate(
        dont_add_these_options,
        prefix_with_dash + "ksp_max_it",
        stringify(params.get<unsigned int>("l_max_its")));
 
    Moose::PetscSupport::setSinglePetscOptionIfAppropriate(dont_add_these_options,
                                                           prefix_with_dash + "ksp_rtol",
                                                           stringify(params.get<Real>("l_tol")));
 
    Moose::PetscSupport::setSinglePetscOptionIfAppropriate(
        dont_add_these_options,
        prefix_with_dash + "ksp_atol",
        stringify(params.get<Real>("l_abs_tol")));
  }
  else
  { // linear eigenvalue problem
    // linear solver
    Moose::PetscSupport::setSinglePetscOptionIfAppropriate(
        dont_add_these_options,
        prefix_with_dash + "st_ksp_max_it",
        stringify(params.get<unsigned int>("l_max_its")));
 
    Moose::PetscSupport::setSinglePetscOptionIfAppropriate(dont_add_these_options,
                                                           prefix_with_dash + "st_ksp_rtol",
                                                           stringify(params.get<Real>("l_tol")));
 
    Moose::PetscSupport::setSinglePetscOptionIfAppropriate(
        dont_add_these_options,
        prefix_with_dash + "st_ksp_atol",
        stringify(params.get<Real>("l_abs_tol")));
  }
}

Referenced by slepcSetOptions().

setWhichEigenPairsOptions()

void Moose::SlepcSupport::setWhichEigenPairsOptions	(	SolverParams &	solver_params,
		const MultiMooseEnum &	dont_add_these_options
	)

Definition at line 339 of file SlepcSupport.C.

{
  switch (solver_params._which_eigen_pairs)
  {
    case Moose::WEP_LARGEST_MAGNITUDE:
      Moose::PetscSupport::setSinglePetscOptionIfAppropriate(dont_add_these_options,
                                                             "-eps_largest_magnitude");
      break;
 
    case Moose::WEP_SMALLEST_MAGNITUDE:
      Moose::PetscSupport::setSinglePetscOptionIfAppropriate(dont_add_these_options,
                                                             "-eps_smallest_magnitude");
      break;
 
    case Moose::WEP_LARGEST_REAL:
      Moose::PetscSupport::setSinglePetscOptionIfAppropriate(dont_add_these_options,
                                                             "-eps_largest_real");
      break;
 
    case Moose::WEP_SMALLEST_REAL:
      Moose::PetscSupport::setSinglePetscOptionIfAppropriate(dont_add_these_options,
                                                             "-eps_smallest_real");
      break;
 
    case Moose::WEP_LARGEST_IMAGINARY:
      Moose::PetscSupport::setSinglePetscOptionIfAppropriate(dont_add_these_options,
                                                             "-eps_largest_imaginary");
      break;
 
    case Moose::WEP_SMALLEST_IMAGINARY:
      Moose::PetscSupport::setSinglePetscOptionIfAppropriate(dont_add_these_options,
                                                             "-eps_smallest_imaginary");
      break;
 
    case Moose::WEP_TARGET_MAGNITUDE:
      Moose::PetscSupport::setSinglePetscOptionIfAppropriate(dont_add_these_options,
                                                             "-eps_target_magnitude");
      break;
 
    case Moose::WEP_TARGET_REAL:
      Moose::PetscSupport::setSinglePetscOptionIfAppropriate(dont_add_these_options,
                                                             "-eps_target_real");
      break;
 
    case Moose::WEP_TARGET_IMAGINARY:
      Moose::PetscSupport::setSinglePetscOptionIfAppropriate(dont_add_these_options,
                                                             "-eps_target_imaginary");
      break;
 
    case Moose::WEP_ALL_EIGENVALUES:
      Moose::PetscSupport::setSinglePetscOptionIfAppropriate(dont_add_these_options, "-eps_all");
      break;
 
    case Moose::WEP_SLEPC_DEFAULT:
      break;
 
    default:
      mooseError("Unknown type of WhichEigenPairs \n");
  }
}

Referenced by slepcSetOptions().

slepcSetOptions()

void Moose::SlepcSupport::slepcSetOptions	(	EigenProblem &	eigen_problem,
		SolverParams &	solver_params,
		const InputParameters &	params
	)

Push all SLEPc/PETSc options into SLEPc/PETSc side.

Options could come from commandline, SolverParams, params, etc.

Definition at line 548 of file SlepcSupport.C.

{
  const auto & dont_add_these_options = eigen_problem.getPetscOptions().dont_add_these_options;
 
  Moose::PetscSupport::petscSetOptions(
      eigen_problem.getPetscOptions(), solver_params, &eigen_problem);
  // Call "SolverTolerances" first, so some solver specific tolerance such as "eps_max_it"
  // can be overriden
  setSlepcEigenSolverTolerances(eigen_problem, solver_params, params);
  setEigenSolverOptions(solver_params, params);
  // when Bx norm postprocessor is provided, we switch off the sign normalization
  if (eigen_problem.bxNormProvided())
    Moose::PetscSupport::setSinglePetscOptionIfAppropriate(
        dont_add_these_options, "-eps_power_sign_normalization", "0", &eigen_problem);
  setEigenProblemOptions(solver_params, eigen_problem.getPetscOptions().dont_add_these_options);
  setWhichEigenPairsOptions(solver_params, eigen_problem.getPetscOptions().dont_add_these_options);
  Moose::PetscSupport::addPetscOptionsFromCommandline();
}

Referenced by EigenProblemSolve::initialSetup().

storeSolveType()

void Moose::SlepcSupport::storeSolveType	(	FEProblemBase &	fe_problem,
		const InputParameters &	params
	)

Set solve type into eigen problem (solverParams)

Definition at line 284 of file SlepcSupport.C.

{
  if (!(dynamic_cast<EigenProblem *>(&fe_problem)))
    return;
 
  if (params.isParamValid("solve_type"))
    for (const auto i : make_range(fe_problem.numNonlinearSystems()))
      fe_problem.solverParams(i)._eigen_solve_type =
          Moose::stringToEnum<Moose::EigenSolveType>(params.get<MooseEnum>("solve_type"));
}

Referenced by EigenProblemSolve::EigenProblemSolve().

Variable Documentation

subspace_factor

const int Moose::SlepcSupport::subspace_factor = 2

Definition at line 32 of file SlepcSupport.C.

Referenced by setEigenProblemSolverParams().