https://mooseframework.inl.gov
SlepcSupport.C
Go to the documentation of this file.
1 //* This file is part of the MOOSE framework
2 //* https://mooseframework.inl.gov
3 //*
4 //* All rights reserved, see COPYRIGHT for full restrictions
5 //* https://github.com/idaholab/moose/blob/master/COPYRIGHT
6 //*
7 //* Licensed under LGPL 2.1, please see LICENSE for details
8 //* https://www.gnu.org/licenses/lgpl-2.1.html
9 
10 #include "libmesh/libmesh_config.h"
11 
12 #ifdef LIBMESH_HAVE_SLEPC
13 
14 #include "SlepcSupport.h"
15 // MOOSE includes
16 #include "InputParameters.h"
17 #include "Conversion.h"
18 #include "EigenProblem.h"
19 #include "FEProblemBase.h"
20 #include "NonlinearEigenSystem.h"
21 #include "libmesh/petsc_vector.h"
22 #include "libmesh/petsc_matrix.h"
23 #include "libmesh/slepc_macro.h"
24 #include "libmesh/auto_ptr.h"
25 #include "petscsnes.h"
26 #include "slepceps.h"
27 
28 using namespace libMesh;
29 
30 namespace Moose
31 {
32 namespace SlepcSupport
33 {
34 
35 const int subspace_factor = 2;
36 
39 {
40  MooseEnum solve_type("POWER ARNOLDI KRYLOVSCHUR JACOBI_DAVIDSON "
41  "NONLINEAR_POWER NEWTON PJFNK PJFNKMO JFNK",
42  "PJFNK");
43  params.set<MooseEnum>("solve_type") = solve_type;
44 
45  params.setDocString("solve_type",
46  "POWER: Power / Inverse / RQI "
47  "ARNOLDI: Arnoldi "
48  "KRYLOVSCHUR: Krylov-Schur "
49  "JACOBI_DAVIDSON: Jacobi-Davidson "
50  "NONLINEAR_POWER: Nonlinear Power "
51  "NEWTON: Newton "
52  "PJFNK: Preconditioned Jacobian-free Newton-Kyrlov"
53  "JFNK: Jacobian-free Newton-Kyrlov"
54  "PJFNKMO: Preconditioned Jacobian-free Newton-Kyrlov with Matrix Only");
55 
56  // When the eigenvalue problems is reformed as a coupled nonlinear system,
57  // we use part of Jacobian as the preconditioning matrix.
58  // Because the difference between the Jacobian and the preconditioning matrix is not small,
59  // the linear solver KSP can not reduce the residual much. After several tests,
60  // we find 1e-2 is a reasonable choice.
61  params.set<Real>("l_tol") = 1e-2;
62 
63  return params;
64 }
65 
68 {
70 
71  // We are solving a Non-Hermitian eigenvalue problem by default
72  MooseEnum eigen_problem_type("HERMITIAN NON_HERMITIAN GEN_HERMITIAN GEN_NON_HERMITIAN "
73  "GEN_INDEFINITE POS_GEN_NON_HERMITIAN SLEPC_DEFAULT",
74  "GEN_NON_HERMITIAN");
75  params.addParam<MooseEnum>(
76  "eigen_problem_type",
77  eigen_problem_type,
78  "Type of the eigenvalue problem we are solving "
79  "HERMITIAN: Hermitian "
80  "NON_HERMITIAN: Non-Hermitian "
81  "GEN_HERMITIAN: Generalized Hermitian "
82  "GEN_NON_HERMITIAN: Generalized Non-Hermitian "
83  "GEN_INDEFINITE: Generalized indefinite Hermitian "
84  "POS_GEN_NON_HERMITIAN: Generalized Non-Hermitian with positive (semi-)definite B "
85  "SLEPC_DEFAULT: Use whatever SLEPC has by default ");
86 
87  // Which eigenvalues are we interested in
88  MooseEnum which_eigen_pairs("LARGEST_MAGNITUDE SMALLEST_MAGNITUDE LARGEST_REAL SMALLEST_REAL "
89  "LARGEST_IMAGINARY SMALLEST_IMAGINARY TARGET_MAGNITUDE TARGET_REAL "
90  "TARGET_IMAGINARY ALL_EIGENVALUES SLEPC_DEFAULT");
91  params.addParam<MooseEnum>("which_eigen_pairs",
92  which_eigen_pairs,
93  "Which eigenvalue pairs to obtain from the solution "
94  "LARGEST_MAGNITUDE "
95  "SMALLEST_MAGNITUDE "
96  "LARGEST_REAL "
97  "SMALLEST_REAL "
98  "LARGEST_IMAGINARY "
99  "SMALLEST_IMAGINARY "
100  "TARGET_MAGNITUDE "
101  "TARGET_REAL "
102  "TARGET_IMAGINARY "
103  "ALL_EIGENVALUES "
104  "SLEPC_DEFAULT ");
105 
106  params.addParam<unsigned int>("n_eigen_pairs", 1, "The number of eigen pairs");
107  params.addParam<unsigned int>("n_basis_vectors", 3, "The dimension of eigen subspaces");
108 
109  params.addParam<Real>("eigen_tol", 1.0e-4, "Relative Tolerance for Eigen Solver");
110  params.addParam<unsigned int>("eigen_max_its", 10000, "Max Iterations for Eigen Solver");
111 
112  params.addParam<Real>("l_abs_tol", 1e-50, "Absolute Tolerances for Linear Solver");
113 
114  params.addParam<unsigned int>("free_power_iterations", 4, "The number of free power iterations");
115 
116  params.addParam<unsigned int>(
117  "extra_power_iterations", 0, "The number of extra free power iterations");
118 
119  params.addParamNamesToGroup(
120  "eigen_problem_type which_eigen_pairs n_eigen_pairs n_basis_vectors eigen_tol eigen_max_its "
121  "free_power_iterations extra_power_iterations",
122  "Eigen Solver");
123  params.addParamNamesToGroup("l_abs_tol", "Linear solver");
124 
125  return params;
126 }
127 
128 void
130  const SolverParams & solver_params,
131  const InputParameters & params)
132 {
133  const auto & dont_add_these_options = eigen_problem.getPetscOptions().dont_add_these_options;
134  const auto prefix_with_dash = '-' + solver_params._prefix;
135 
136  mooseAssert(solver_params._solver_sys_num != libMesh::invalid_uint,
137  "The solver system number must be initialized");
138 
140  prefix_with_dash + "eps_tol",
141  stringify(params.get<Real>("eigen_tol")));
142 
144  dont_add_these_options,
145  prefix_with_dash + "eps_max_it",
146  stringify(params.get<unsigned int>("eigen_max_its")));
147 
148  // if it is a nonlinear eigenvalue solver, we need to set tolerances for nonlinear solver and
149  // linear solver
150  if (eigen_problem.isNonlinearEigenvalueSolver(solver_params._solver_sys_num))
151  {
152  // nonlinear solver tolerances
154  dont_add_these_options,
155  prefix_with_dash + "snes_max_it",
156  stringify(params.get<unsigned int>("nl_max_its")));
157 
159  dont_add_these_options,
160  prefix_with_dash + "snes_max_funcs",
161  stringify(params.get<unsigned int>("nl_max_funcs")));
162 
164  dont_add_these_options,
165  prefix_with_dash + "snes_atol",
166  stringify(params.get<Real>("nl_abs_tol")));
167 
169  dont_add_these_options,
170  prefix_with_dash + "snes_rtol",
171  stringify(params.get<Real>("nl_rel_tol")));
172 
174  dont_add_these_options,
175  prefix_with_dash + "snes_stol",
176  stringify(params.get<Real>("nl_rel_step_tol")));
177 
178  // linear solver
180  dont_add_these_options,
181  prefix_with_dash + "ksp_max_it",
182  stringify(params.get<unsigned int>("l_max_its")));
183 
185  prefix_with_dash + "ksp_rtol",
186  stringify(params.get<Real>("l_tol")));
187 
189  dont_add_these_options,
190  prefix_with_dash + "ksp_atol",
191  stringify(params.get<Real>("l_abs_tol")));
192  }
193  else
194  { // linear eigenvalue problem
195  // linear solver
197  dont_add_these_options,
198  prefix_with_dash + "st_ksp_max_it",
199  stringify(params.get<unsigned int>("l_max_its")));
200 
202  prefix_with_dash + "st_ksp_rtol",
203  stringify(params.get<Real>("l_tol")));
204 
206  dont_add_these_options,
207  prefix_with_dash + "st_ksp_atol",
208  stringify(params.get<Real>("l_abs_tol")));
209  }
210 }
211 
212 void
214 {
215  for (const auto i : make_range(eigen_problem.numNonlinearSystems()))
216  {
217  const std::string & eigen_problem_type = params.get<MooseEnum>("eigen_problem_type");
218  if (!eigen_problem_type.empty())
219  eigen_problem.solverParams(i)._eigen_problem_type =
220  Moose::stringToEnum<Moose::EigenProblemType>(eigen_problem_type);
221  else
222  mooseError("Have to specify a valid eigen problem type");
223 
224  const std::string & which_eigen_pairs = params.get<MooseEnum>("which_eigen_pairs");
225  if (!which_eigen_pairs.empty())
226  eigen_problem.solverParams(i)._which_eigen_pairs =
227  Moose::stringToEnum<Moose::WhichEigenPairs>(which_eigen_pairs);
228 
229  // Set necessary parameters used in EigenSystem::solve(),
230  // i.e. the number of requested eigenpairs nev and the number
231  // of basis vectors ncv used in the solution algorithm. Note that
232  // ncv >= nev must hold and ncv >= 2*nev is recommended
233  unsigned int n_eigen_pairs = params.get<unsigned int>("n_eigen_pairs");
234  unsigned int n_basis_vectors = params.get<unsigned int>("n_basis_vectors");
235 
236  eigen_problem.setNEigenPairsRequired(n_eigen_pairs);
237 
238  eigen_problem.es().parameters.set<unsigned int>("eigenpairs") = n_eigen_pairs;
239 
240  // If the subspace dimension is too small, we increase it automatically
241  if (subspace_factor * n_eigen_pairs > n_basis_vectors)
242  {
243  n_basis_vectors = subspace_factor * n_eigen_pairs;
244  mooseWarning(
245  "Number of subspaces in Eigensolver is changed by moose because the value you set "
246  "is too small");
247  }
248 
249  eigen_problem.es().parameters.set<unsigned int>("basis vectors") = n_basis_vectors;
250 
251  // Operators A and B are formed as shell matrices
252  eigen_problem.solverParams(i)._eigen_matrix_free = params.get<bool>("matrix_free");
253 
254  // Preconditioning is formed as a shell matrix
255  eigen_problem.solverParams(i)._precond_matrix_free = params.get<bool>("precond_matrix_free");
256 
257  if (params.get<MooseEnum>("solve_type") == "PJFNK")
258  {
259  eigen_problem.solverParams(i)._eigen_matrix_free = true;
260  }
261  if (params.get<MooseEnum>("solve_type") == "JFNK")
262  {
263  eigen_problem.solverParams(i)._eigen_matrix_free = true;
264  eigen_problem.solverParams(i)._precond_matrix_free = true;
265  }
266  // We need matrices so that we can implement residual evaluations
267  if (params.get<MooseEnum>("solve_type") == "PJFNKMO")
268  {
269  eigen_problem.solverParams(i)._eigen_matrix_free = true;
270  eigen_problem.solverParams(i)._precond_matrix_free = false;
271  eigen_problem.solverParams(i)._eigen_matrix_vector_mult = true;
272  // By default, we need to form full matrices, otherwise residual
273  // evaluations will not be accurate
274  eigen_problem.setCoupling(Moose::COUPLING_FULL);
275  }
276  }
277 
278  eigen_problem.constantMatrices(params.get<bool>("constant_matrices"));
279 
280  if (eigen_problem.constantMatrices() && params.get<MooseEnum>("solve_type") != "PJFNKMO")
281  {
282  mooseError("constant_matrices flag is only valid for solve type: PJFNKMO");
283  }
284 }
285 
286 void
287 storeSolveType(FEProblemBase & fe_problem, const InputParameters & params)
288 {
289  if (!(dynamic_cast<EigenProblem *>(&fe_problem)))
290  return;
291 
292  if (params.isParamValid("solve_type"))
293  for (const auto i : make_range(fe_problem.numNonlinearSystems()))
294  fe_problem.solverParams(i)._eigen_solve_type =
295  Moose::stringToEnum<Moose::EigenSolveType>(params.get<MooseEnum>("solve_type"));
296 }
297 
298 void
299 setEigenProblemOptions(SolverParams & solver_params, const MultiMooseEnum & dont_add_these_options)
300 {
301  switch (solver_params._eigen_problem_type)
302  {
305  "-eps_hermitian");
306  break;
307 
310  "-eps_non_hermitian");
311  break;
312 
315  "-eps_gen_hermitian");
316  break;
317 
320  "-eps_gen_indefinite");
321  break;
322 
325  "-eps_gen_non_hermitian");
326  break;
327 
330  "-eps_pos_gen_non_hermitian");
331  break;
332 
334  break;
335 
336  default:
337  mooseError("Unknown eigen solver type \n");
338  }
339 }
340 
341 void
343  const MultiMooseEnum & dont_add_these_options)
344 {
345  switch (solver_params._which_eigen_pairs)
346  {
349  "-eps_largest_magnitude");
350  break;
351 
354  "-eps_smallest_magnitude");
355  break;
356 
359  "-eps_largest_real");
360  break;
361 
364  "-eps_smallest_real");
365  break;
366 
369  "-eps_largest_imaginary");
370  break;
371 
374  "-eps_smallest_imaginary");
375  break;
376 
379  "-eps_target_magnitude");
380  break;
381 
384  "-eps_target_real");
385  break;
386 
389  "-eps_target_imaginary");
390  break;
391 
393  Moose::PetscSupport::setSinglePetscOptionIfAppropriate(dont_add_these_options, "-eps_all");
394  break;
395 
397  break;
398 
399  default:
400  mooseError("Unknown type of WhichEigenPairs \n");
401  }
402 }
403 
404 void
405 setFreeNonlinearPowerIterations(unsigned int free_power_iterations)
406 {
407  Moose::PetscSupport::setSinglePetscOption("-eps_power_update", "0");
408  Moose::PetscSupport::setSinglePetscOption("-snes_max_it", "2");
409  // During each power iteration, we want solver converged unless linear solver does not
410  // work. We here use a really loose tolerance for this purpose.
411  // -snes_no_convergence_test is a perfect option, but it was removed from PETSc
412  Moose::PetscSupport::setSinglePetscOption("-snes_rtol", "0.99999999999");
413  Moose::PetscSupport::setSinglePetscOption("-eps_max_it", stringify(free_power_iterations));
414  // We always want the number of free power iterations respected so we don't want to stop early if
415  // we've satisfied a convergence criterion. Consequently we make this tolerance very tight
416  Moose::PetscSupport::setSinglePetscOption("-eps_tol", "1e-50");
417 }
418 
419 void
421 {
422  Moose::PetscSupport::setSinglePetscOption("-eps_power_update", "1");
423  Moose::PetscSupport::setSinglePetscOption("-eps_max_it", "1");
425  stringify(params.get<unsigned int>("nl_max_its")));
427  stringify(params.get<Real>("nl_rel_tol")));
428  Moose::PetscSupport::setSinglePetscOption("-eps_tol", stringify(params.get<Real>("eigen_tol")));
429 }
430 
431 void
432 setNewtonPetscOptions(SolverParams & solver_params, const InputParameters & params)
433 {
434 #if !SLEPC_VERSION_LESS_THAN(3, 8, 0) || !PETSC_VERSION_RELEASE
435  // Whether or not we need to involve an initial inverse power
436  bool initial_power = params.get<bool>("_newton_inverse_power");
437 
438  Moose::PetscSupport::setSinglePetscOption("-eps_type", "power");
439  Moose::PetscSupport::setSinglePetscOption("-eps_power_nonlinear", "1");
440  Moose::PetscSupport::setSinglePetscOption("-eps_power_update", "1");
441  // Only one outer iteration in EPS is allowed when Newton/PJFNK/JFNK
442  // is used as the eigen solver
443  Moose::PetscSupport::setSinglePetscOption("-eps_max_it", "1");
444  if (initial_power)
445  {
446  Moose::PetscSupport::setSinglePetscOption("-init_eps_power_snes_max_it", "1");
447  Moose::PetscSupport::setSinglePetscOption("-init_eps_power_ksp_rtol", "1e-2");
449  "-init_eps_max_it", stringify(params.get<unsigned int>("free_power_iterations")));
450  }
451  Moose::PetscSupport::setSinglePetscOption("-eps_target_magnitude", "");
452  if (solver_params._eigen_matrix_free)
453  {
454  Moose::PetscSupport::setSinglePetscOption("-snes_mf_operator", "1");
455  if (initial_power)
456  Moose::PetscSupport::setSinglePetscOption("-init_eps_power_snes_mf_operator", "1");
457  }
458  else
459  {
460  Moose::PetscSupport::setSinglePetscOption("-snes_mf_operator", "0");
461  if (initial_power)
462  Moose::PetscSupport::setSinglePetscOption("-init_eps_power_snes_mf_operator", "0");
463  }
464 #if PETSC_RELEASE_LESS_THAN(3, 13, 0)
465  Moose::PetscSupport::setSinglePetscOption("-st_type", "sinvert");
466  if (initial_power)
467  Moose::PetscSupport::setSinglePetscOption("-init_st_type", "sinvert");
468 #endif
469 #else
470  mooseError("Newton-based eigenvalue solver requires SLEPc 3.7.3 or higher");
471 #endif
472 }
473 
474 void
476 {
477 #if !SLEPC_VERSION_LESS_THAN(3, 8, 0) || !PETSC_VERSION_RELEASE
478  Moose::PetscSupport::setSinglePetscOption("-eps_type", "power");
479  Moose::PetscSupport::setSinglePetscOption("-eps_power_nonlinear", "1");
480  Moose::PetscSupport::setSinglePetscOption("-eps_target_magnitude", "");
481  if (solver_params._eigen_matrix_free)
482  Moose::PetscSupport::setSinglePetscOption("-snes_mf_operator", "1");
483  else
484  Moose::PetscSupport::setSinglePetscOption("-snes_mf_operator", "0");
485 
486 #if PETSC_RELEASE_LESS_THAN(3, 13, 0)
487  Moose::PetscSupport::setSinglePetscOption("-st_type", "sinvert");
488 #endif
489 #else
490  mooseError("Nonlinear Inverse Power requires SLEPc 3.7.3 or higher");
491 #endif
492 }
493 
494 void
495 setEigenSolverOptions(SolverParams & solver_params, const InputParameters & params)
496 {
497  // Avoid unused variable warnings when you have SLEPc but not PETSc-dev.
498  libmesh_ignore(params);
499 
500  switch (solver_params._eigen_solve_type)
501  {
502  case Moose::EST_POWER:
503  Moose::PetscSupport::setSinglePetscOption("-eps_type", "power");
504  break;
505 
506  case Moose::EST_ARNOLDI:
507  Moose::PetscSupport::setSinglePetscOption("-eps_type", "arnoldi");
508  break;
509 
511  Moose::PetscSupport::setSinglePetscOption("-eps_type", "krylovschur");
512  break;
513 
516  break;
517 
519  setNonlinearPowerOptions(solver_params);
520  break;
521 
522  case Moose::EST_NEWTON:
523  setNewtonPetscOptions(solver_params, params);
524  break;
525 
526  case Moose::EST_PJFNK:
527  solver_params._eigen_matrix_free = true;
528  solver_params._customized_pc_for_eigen = false;
529  setNewtonPetscOptions(solver_params, params);
530  break;
531 
532  case Moose::EST_JFNK:
533  solver_params._eigen_matrix_free = true;
534  solver_params._customized_pc_for_eigen = true;
535  setNewtonPetscOptions(solver_params, params);
536  break;
537 
538  case Moose::EST_PJFNKMO:
539  solver_params._eigen_matrix_free = true;
540  solver_params._customized_pc_for_eigen = false;
541  solver_params._eigen_matrix_vector_mult = true;
542  setNewtonPetscOptions(solver_params, params);
543  break;
544 
545  default:
546  mooseError("Unknown eigen solver type \n");
547  }
548 }
549 
550 void
552  SolverParams & solver_params,
553  const InputParameters & params)
554 {
555  const auto & dont_add_these_options = eigen_problem.getPetscOptions().dont_add_these_options;
556 
558  eigen_problem.getPetscOptions(), solver_params, &eigen_problem);
559  // Call "SolverTolerances" first, so some solver specific tolerance such as "eps_max_it"
560  // can be overriden
561  setSlepcEigenSolverTolerances(eigen_problem, solver_params, params);
562  setEigenSolverOptions(solver_params, params);
563  // when Bx norm postprocessor is provided, we switch off the sign normalization
564  if (eigen_problem.bxNormProvided())
566  dont_add_these_options, "-eps_power_sign_normalization", "0", &eigen_problem);
567  setEigenProblemOptions(solver_params, eigen_problem.getPetscOptions().dont_add_these_options);
568  setWhichEigenPairsOptions(solver_params, eigen_problem.getPetscOptions().dont_add_these_options);
570 }
571 
572 // For matrices A and B
573 PetscErrorCode
574 mooseEPSFormMatrices(EigenProblem & eigen_problem, EPS eps, Vec x, void * ctx)
575 {
576  ST st;
577  Mat A, B;
578  PetscBool aisshell, bisshell;
580 
581  if (eigen_problem.constantMatrices() && eigen_problem.wereMatricesFormed())
582  PetscFunctionReturn(PETSC_SUCCESS);
583 
584  if (eigen_problem.onLinearSolver())
585  // We reach here during linear iteration when solve type is PJFNKMO.
586  // We will use the matrices assembled at the beginning of this Newton
587  // iteration for the following residual evaluation.
588  PetscFunctionReturn(PETSC_SUCCESS);
589 
590  NonlinearEigenSystem & eigen_nl = eigen_problem.getCurrentNonlinearEigenSystem();
591  auto & sys = eigen_nl.sys();
592  SNES snes = eigen_nl.getSNES();
593  // Rest ST state so that we can retrieve matrices
594  LibmeshPetscCallQ(EPSGetST(eps, &st));
595  LibmeshPetscCallQ(STResetMatrixState(st));
596  LibmeshPetscCallQ(EPSGetOperators(eps, &A, &B));
597  LibmeshPetscCallQ(PetscObjectTypeCompare((PetscObject)A, MATSHELL, &aisshell));
598  LibmeshPetscCallQ(PetscObjectTypeCompare((PetscObject)B, MATSHELL, &bisshell));
599  if (aisshell || bisshell)
600  {
601  SETERRQ(PetscObjectComm((PetscObject)eps),
602  PETSC_ERR_ARG_INCOMP,
603  "A and B matrices can not be shell matrices when using PJFNKMO \n");
604  }
605  // Form A and B
606  std::vector<Mat> mats = {A, B};
607  std::vector<SparseMatrix<Number> *> libmesh_mats = {&sys.get_matrix_A(), &sys.get_matrix_B()};
609  snes, x, mats, libmesh_mats, ctx, {eigen_nl.nonEigenMatrixTag(), eigen_nl.eigenMatrixTag()});
610  eigen_problem.wereMatricesFormed(true);
611  PetscFunctionReturn(PETSC_SUCCESS);
612 }
613 
614 namespace
615 {
616 void
617 updateCurrentLocalSolution(CondensedEigenSystem & sys, Vec x)
618 {
619  auto & dof_map = sys.get_dof_map();
620 
621  PetscVector<Number> X_global(x, sys.comm());
622 
623  if (dof_map.n_constrained_dofs())
624  {
625  sys.copy_sub_to_super(X_global, *sys.solution);
626  // Set the constrained dof values
627  dof_map.enforce_constraints_exactly(sys);
628  sys.update();
629  }
630  else
631  {
632  PetscVector<Number> & X_sys = *cast_ptr<PetscVector<Number> *>(sys.solution.get());
633 
634  // Use the system's update() to get a good local version of the
635  // parallel solution. This operation does not modify the incoming
636  // "x" vector, it only localizes information from "x" into
637  // sys.current_local_solution.
638  X_global.swap(X_sys);
639  sys.update();
640  X_global.swap(X_sys);
641  }
642 }
643 
644 std::unique_ptr<NumericVector<Number>>
645 createWrappedResidual(CondensedEigenSystem & sys, Vec r)
646 {
647  auto & dof_map = sys.get_dof_map();
648 
649  if (dof_map.n_constrained_dofs())
650  return sys.solution->zero_clone();
651  else
652  {
653  auto R = std::make_unique<PetscVector<Number>>(r, sys.comm());
654  R->zero();
655  return R;
656  }
657 }
658 
659 void
660 evaluateResidual(EigenProblem & eigen_problem, Vec x, Vec r, TagID tag)
661 {
662  auto & nl = eigen_problem.getCurrentNonlinearEigenSystem();
663  auto & sys = nl.sys();
664  auto & dof_map = sys.get_dof_map();
665 
666  updateCurrentLocalSolution(sys, x);
667  auto R = createWrappedResidual(sys, r);
668 
669  eigen_problem.computeResidualTag(*sys.current_local_solution.get(), *R, tag);
670 
671  R->close();
672 
673  if (dof_map.n_constrained_dofs())
674  {
675  PetscVector<Number> sub_r(r, sys.comm());
676  sys.copy_super_to_sub(*R, sub_r);
677  }
678 }
679 }
680 
681 void
683  SNES /*snes*/, Vec x, Mat eigen_mat, SparseMatrix<Number> & all_dofs_mat, void * ctx, TagID tag)
684 {
685  EigenProblem * eigen_problem = static_cast<EigenProblem *>(ctx);
686  auto & nl = eigen_problem->getCurrentNonlinearEigenSystem();
687  auto & sys = nl.sys();
688  auto & dof_map = sys.get_dof_map();
689 
690 #ifndef NDEBUG
691  auto & petsc_all_dofs_mat = cast_ref<PetscMatrix<Number> &>(all_dofs_mat);
692  mooseAssert(
693  !dof_map.n_constrained_dofs() == (eigen_mat == petsc_all_dofs_mat.mat()),
694  "If we do not have constrained dofs, then eigen_mat and all_dofs_mat should be the same. "
695  "Conversely, if we do have constrained dofs, they must be different");
696 #endif
697 
698  updateCurrentLocalSolution(sys, x);
699 
700  if (!eigen_problem->constJacobian())
701  all_dofs_mat.zero();
702 
703  eigen_problem->computeJacobianTag(*sys.current_local_solution.get(), all_dofs_mat, tag);
704 
705  if (dof_map.n_constrained_dofs())
706  {
707  PetscMatrix<Number> wrapped_eigen_mat(eigen_mat, sys.comm());
708  sys.copy_super_to_sub(all_dofs_mat, wrapped_eigen_mat);
709  }
710 }
711 
712 void
714  Vec x,
715  std::vector<Mat> & eigen_mats,
716  std::vector<SparseMatrix<Number> *> & all_dofs_mats,
717  void * ctx,
718  const std::set<TagID> & tags)
719 {
720  EigenProblem * eigen_problem = static_cast<EigenProblem *>(ctx);
721  auto & nl = eigen_problem->getCurrentNonlinearEigenSystem();
722  auto & sys = nl.sys();
723  auto & dof_map = sys.get_dof_map();
724 
725 #ifndef NDEBUG
726  for (const auto i : index_range(eigen_mats))
727  mooseAssert(!dof_map.n_constrained_dofs() ==
728  (eigen_mats[i] == cast_ptr<PetscMatrix<Number> *>(all_dofs_mats[i])->mat()),
729  "If we do not have constrained dofs, then mat and libmesh_mat should be the same. "
730  "Conversely, if we do have constrained dofs, they must be different");
731 #endif
732 
733  updateCurrentLocalSolution(sys, x);
734 
735  for (auto * const all_dofs_mat : all_dofs_mats)
736  if (!eigen_problem->constJacobian())
737  all_dofs_mat->zero();
738 
739  eigen_problem->computeMatricesTags(*sys.current_local_solution.get(), all_dofs_mats, tags);
740 
741  if (dof_map.n_constrained_dofs())
742  for (const auto i : index_range(eigen_mats))
743  {
744  PetscMatrix<Number> wrapped_eigen_mat(eigen_mats[i], sys.comm());
745  sys.copy_super_to_sub(*all_dofs_mats[i], wrapped_eigen_mat);
746  }
747 }
748 
749 PetscErrorCode
751 {
752  PetscErrorCode (*func)(SNES, Vec, Vec, void *);
753  void * fctx;
755 
756  EigenProblem * eigen_problem = static_cast<EigenProblem *>(ctx);
757  NonlinearEigenSystem & eigen_nl = eigen_problem->getCurrentNonlinearEigenSystem();
758  SNES snes = eigen_nl.getSNES();
759 
760  eigen_problem->onLinearSolver(true);
761 
762  LibmeshPetscCallQ(SNESGetFunction(snes, NULL, &func, &fctx));
763  if (fctx != ctx)
764  {
765  SETERRQ(
766  PetscObjectComm((PetscObject)snes), PETSC_ERR_ARG_INCOMP, "Contexts are not consistent \n");
767  }
768  LibmeshPetscCallQ((*func)(snes, x, r, ctx));
769 
770  eigen_problem->onLinearSolver(false);
771 
772  PetscFunctionReturn(PETSC_SUCCESS);
773 }
774 
775 PetscErrorCode
776 mooseSlepcEigenFormJacobianA(SNES snes, Vec x, Mat jac, Mat pc, void * ctx)
777 {
778  PetscBool jisshell, pisshell;
779  PetscBool jismffd;
780 
782 
783  EigenProblem * eigen_problem = static_cast<EigenProblem *>(ctx);
784  NonlinearEigenSystem & eigen_nl = eigen_problem->getCurrentNonlinearEigenSystem();
785  auto & sys = eigen_nl.sys();
786 
787  // If both jacobian and preconditioning are shell matrices,
788  // and then assemble them and return
789  LibmeshPetscCallQ(PetscObjectTypeCompare((PetscObject)jac, MATSHELL, &jisshell));
790  LibmeshPetscCallQ(PetscObjectTypeCompare((PetscObject)jac, MATMFFD, &jismffd));
791 
792  if (jismffd && eigen_problem->solverParams(eigen_nl.number())._eigen_matrix_vector_mult)
793  {
796 
797  EPS eps = eigen_nl.getEPS();
798 
799  LibmeshPetscCallQ(mooseEPSFormMatrices(*eigen_problem, eps, x, ctx));
800 
801  if (pc != jac)
802  {
803  LibmeshPetscCallQ(MatAssemblyBegin(jac, MAT_FINAL_ASSEMBLY));
804  LibmeshPetscCallQ(MatAssemblyEnd(jac, MAT_FINAL_ASSEMBLY));
805  }
806  PetscFunctionReturn(PETSC_SUCCESS);
807  }
808 
809  LibmeshPetscCallQ(PetscObjectTypeCompare((PetscObject)pc, MATSHELL, &pisshell));
810  if ((jisshell || jismffd) && pisshell)
811  {
812  // Just assemble matrices and return
813  LibmeshPetscCallQ(MatAssemblyBegin(jac, MAT_FINAL_ASSEMBLY));
814  LibmeshPetscCallQ(MatAssemblyBegin(pc, MAT_FINAL_ASSEMBLY));
815  LibmeshPetscCallQ(MatAssemblyEnd(jac, MAT_FINAL_ASSEMBLY));
816  LibmeshPetscCallQ(MatAssemblyEnd(pc, MAT_FINAL_ASSEMBLY));
817 
818  PetscFunctionReturn(PETSC_SUCCESS);
819  }
820 
821  // Jacobian and precond matrix are the same
822  if (jac == pc)
823  {
824  if (!pisshell)
826  snes, x, pc, sys.get_matrix_A(), ctx, eigen_nl.precondMatrixTag());
827 
828  PetscFunctionReturn(PETSC_SUCCESS);
829  }
830  else
831  {
832  if (!jisshell && !jismffd && !pisshell) // We need to form both Jacobian and precond matrix
833  {
834  std::vector<Mat> mats = {jac, pc};
835  std::vector<SparseMatrix<Number> *> libmesh_mats = {&sys.get_matrix_A(),
836  &sys.get_precond_matrix()};
837  std::set<TagID> tags = {eigen_nl.nonEigenMatrixTag(), eigen_nl.precondMatrixTag()};
838  moosePetscSNESFormMatricesTags(snes, x, mats, libmesh_mats, ctx, tags);
839  PetscFunctionReturn(PETSC_SUCCESS);
840  }
841  if (!pisshell) // We need to form only precond matrix
842  {
844  snes, x, pc, sys.get_precond_matrix(), ctx, eigen_nl.precondMatrixTag());
845  LibmeshPetscCallQ(MatAssemblyBegin(jac, MAT_FINAL_ASSEMBLY));
846  LibmeshPetscCallQ(MatAssemblyEnd(jac, MAT_FINAL_ASSEMBLY));
847  PetscFunctionReturn(PETSC_SUCCESS);
848  }
849  if (!jisshell && !jismffd) // We need to form only Jacobian matrix
850  {
852  snes, x, jac, sys.get_matrix_A(), ctx, eigen_nl.nonEigenMatrixTag());
853  LibmeshPetscCallQ(MatAssemblyBegin(pc, MAT_FINAL_ASSEMBLY));
854  LibmeshPetscCallQ(MatAssemblyEnd(pc, MAT_FINAL_ASSEMBLY));
855  PetscFunctionReturn(PETSC_SUCCESS);
856  }
857  }
858  PetscFunctionReturn(PETSC_SUCCESS);
859 }
860 
861 PetscErrorCode
862 mooseSlepcEigenFormJacobianB(SNES snes, Vec x, Mat jac, Mat pc, void * ctx)
863 {
864  PetscBool jshell, pshell;
865  PetscBool jismffd;
866 
868 
869  EigenProblem * eigen_problem = static_cast<EigenProblem *>(ctx);
870  NonlinearEigenSystem & eigen_nl = eigen_problem->getCurrentNonlinearEigenSystem();
871  auto & sys = eigen_nl.sys();
872 
873  // If both jacobian and preconditioning are shell matrices,
874  // and then assemble them and return
875  LibmeshPetscCallQ(PetscObjectTypeCompare((PetscObject)jac, MATSHELL, &jshell));
876  LibmeshPetscCallQ(PetscObjectTypeCompare((PetscObject)jac, MATMFFD, &jismffd));
877  LibmeshPetscCallQ(PetscObjectTypeCompare((PetscObject)pc, MATSHELL, &pshell));
878  if ((jshell || jismffd) && pshell)
879  {
880  // Just assemble matrices and return
881  LibmeshPetscCallQ(MatAssemblyBegin(jac, MAT_FINAL_ASSEMBLY));
882  LibmeshPetscCallQ(MatAssemblyBegin(pc, MAT_FINAL_ASSEMBLY));
883  LibmeshPetscCallQ(MatAssemblyEnd(jac, MAT_FINAL_ASSEMBLY));
884  LibmeshPetscCallQ(MatAssemblyEnd(pc, MAT_FINAL_ASSEMBLY));
885 
886  PetscFunctionReturn(PETSC_SUCCESS);
887  }
888 
889  if (jac != pc && (!jshell && !jshell))
890  SETERRQ(PetscObjectComm((PetscObject)snes),
891  PETSC_ERR_ARG_INCOMP,
892  "Jacobian and precond matrices should be the same for eigen kernels \n");
893 
894  moosePetscSNESFormMatrixTag(snes, x, pc, sys.get_matrix_B(), ctx, eigen_nl.eigenMatrixTag());
895 
896  if (eigen_problem->negativeSignEigenKernel())
897  {
898  LibmeshPetscCallQ(MatScale(pc, -1.));
899  }
900 
901  PetscFunctionReturn(PETSC_SUCCESS);
902 }
903 
904 void
905 moosePetscSNESFormFunction(SNES /*snes*/, Vec x, Vec r, void * ctx, TagID tag)
906 {
907  EigenProblem * eigen_problem = static_cast<EigenProblem *>(ctx);
908  evaluateResidual(*eigen_problem, x, r, tag);
909 }
910 
911 PetscErrorCode
912 mooseSlepcEigenFormFunctionA(SNES snes, Vec x, Vec r, void * ctx)
913 {
915 
916  EigenProblem * eigen_problem = static_cast<EigenProblem *>(ctx);
917  NonlinearEigenSystem & eigen_nl = eigen_problem->getCurrentNonlinearEigenSystem();
918 
919  if (eigen_problem->solverParams(eigen_nl.number())._eigen_matrix_vector_mult &&
920  (eigen_problem->onLinearSolver() || eigen_problem->constantMatrices()))
921  {
922  EPS eps = eigen_nl.getEPS();
923  Mat A;
924  ST st;
925 
926  LibmeshPetscCallQ(mooseEPSFormMatrices(*eigen_problem, eps, x, ctx));
927 
928  // Rest ST state so that we can restrieve matrices
929  LibmeshPetscCallQ(EPSGetST(eps, &st));
930  LibmeshPetscCallQ(STResetMatrixState(st));
931  LibmeshPetscCallQ(EPSGetOperators(eps, &A, NULL));
932 
933  LibmeshPetscCallQ(MatMult(A, x, r));
934 
935  PetscFunctionReturn(PETSC_SUCCESS);
936  }
937 
938  moosePetscSNESFormFunction(snes, x, r, ctx, eigen_nl.nonEigenVectorTag());
939 
940  PetscFunctionReturn(PETSC_SUCCESS);
941 }
942 
943 PetscErrorCode
944 mooseSlepcEigenFormFunctionB(SNES snes, Vec x, Vec r, void * ctx)
945 {
947 
948  EigenProblem * eigen_problem = static_cast<EigenProblem *>(ctx);
949  NonlinearEigenSystem & eigen_nl = eigen_problem->getCurrentNonlinearEigenSystem();
950 
951  if (eigen_problem->solverParams(eigen_nl.number())._eigen_matrix_vector_mult &&
952  (eigen_problem->onLinearSolver() || eigen_problem->constantMatrices()))
953  {
954  EPS eps = eigen_nl.getEPS();
955  ST st;
956  Mat B;
957 
958  LibmeshPetscCallQ(mooseEPSFormMatrices(*eigen_problem, eps, x, ctx));
959 
960  // Rest ST state so that we can restrieve matrices
961  LibmeshPetscCallQ(EPSGetST(eps, &st));
962  LibmeshPetscCallQ(STResetMatrixState(st));
963  LibmeshPetscCallQ(EPSGetOperators(eps, NULL, &B));
964 
965  LibmeshPetscCallQ(MatMult(B, x, r));
966 
967  if (eigen_problem->bxNormProvided())
968  {
969  // User has provided a postprocessor. We need it updated
970  updateCurrentLocalSolution(eigen_nl.sys(), x);
971  eigen_problem->execute(EXEC_LINEAR);
972  }
973  }
974  else
975  moosePetscSNESFormFunction(snes, x, r, ctx, eigen_nl.eigenVectorTag());
976 
977  if (eigen_problem->negativeSignEigenKernel())
978  {
979  LibmeshPetscCallQ(VecScale(r, -1.));
980  }
981 
982  PetscFunctionReturn(PETSC_SUCCESS);
983 }
984 
985 PetscErrorCode
986 mooseSlepcEigenFormFunctionAB(SNES /*snes*/, Vec x, Vec Ax, Vec Bx, void * ctx)
987 {
989 
990  EigenProblem * eigen_problem = static_cast<EigenProblem *>(ctx);
991  NonlinearEigenSystem & eigen_nl = eigen_problem->getCurrentNonlinearEigenSystem();
992  auto & sys = eigen_nl.sys();
993  auto & dof_map = sys.get_dof_map();
994 
995  if (eigen_problem->solverParams(eigen_nl.number())._eigen_matrix_vector_mult &&
996  (eigen_problem->onLinearSolver() || eigen_problem->constantMatrices()))
997  {
998  EPS eps = eigen_nl.getEPS();
999  ST st;
1000  Mat A, B;
1001 
1002  LibmeshPetscCallQ(mooseEPSFormMatrices(*eigen_problem, eps, x, ctx));
1003 
1004  // Rest ST state so that we can restrieve matrices
1005  LibmeshPetscCallQ(EPSGetST(eps, &st));
1006  LibmeshPetscCallQ(STResetMatrixState(st));
1007 
1008  LibmeshPetscCallQ(EPSGetOperators(eps, &A, &B));
1009 
1010  LibmeshPetscCallQ(MatMult(A, x, Ax));
1011  LibmeshPetscCallQ(MatMult(B, x, Bx));
1012 
1013  if (eigen_problem->negativeSignEigenKernel())
1014  LibmeshPetscCallQ(VecScale(Bx, -1.));
1015 
1016  if (eigen_problem->bxNormProvided())
1017  {
1018  // User has provided a postprocessor. We need it updated
1019  updateCurrentLocalSolution(sys, x);
1020  eigen_problem->execute(EXEC_LINEAR);
1021  }
1022 
1023  PetscFunctionReturn(PETSC_SUCCESS);
1024  }
1025 
1026  updateCurrentLocalSolution(sys, x);
1027  auto AX = createWrappedResidual(sys, Ax);
1028  auto BX = createWrappedResidual(sys, Bx);
1029 
1030  eigen_problem->computeResidualAB(*sys.current_local_solution.get(),
1031  *AX,
1032  *BX,
1033  eigen_nl.nonEigenVectorTag(),
1034  eigen_nl.eigenVectorTag());
1035 
1036  AX->close();
1037  BX->close();
1038 
1039  if (dof_map.n_constrained_dofs())
1040  {
1041  PetscVector<Number> sub_Ax(Ax, sys.comm());
1042  sys.copy_super_to_sub(*AX, sub_Ax);
1043  PetscVector<Number> sub_Bx(Bx, sys.comm());
1044  sys.copy_super_to_sub(*BX, sub_Bx);
1045  }
1046 
1047  if (eigen_problem->negativeSignEigenKernel())
1048  LibmeshPetscCallQ(VecScale(Bx, -1.));
1049 
1050  PetscFunctionReturn(PETSC_SUCCESS);
1051 }
1052 
1053 PetscErrorCode
1054 mooseSlepcEigenFormNorm(SNES /*snes*/, Vec /*Bx*/, PetscReal * norm, void * ctx)
1055 {
1057  auto * const eigen_problem = static_cast<EigenProblem *>(ctx);
1058  *norm = eigen_problem->formNorm();
1059  PetscFunctionReturn(PETSC_SUCCESS);
1060 }
1061 
1062 void
1063 attachCallbacksToMat(EigenProblem & eigen_problem, Mat mat, bool eigen)
1064 {
1065  // Recall that we are solving the potentially nonlinear problem:
1066  // F(x) = A(x) - \lambda B(x) = 0
1067  //
1068  // To solve this, we can use Newton's method: J \Delta x = -F
1069  // Generally we will approximate J using matrix free methods. However, in order to solve the
1070  // linearized system efficiently, we typically will need preconditioning. Typically we will build
1071  // the preconditioner only from A, but we also have the option to include information from B
1072 
1073  // Attach the Jacobian computation function. If \p mat is the "eigen" matrix corresponding to B,
1074  // then attach our JacobianB computation routine, else the matrix corresponds to A, and we attach
1075  // the JacobianA computation routine
1076  LibmeshPetscCallA(
1077  eigen_problem.comm().get(),
1078  PetscObjectComposeFunction((PetscObject)mat,
1079  "formJacobian",
1082 
1083  // Attach the residual computation function. If \p mat is the "eigen" matrix corresponding to B,
1084  // then attach our FunctionB computation routine, else the matrix corresponds to A, and we attach
1085  // the FunctionA computation routine
1086  LibmeshPetscCallA(
1087  eigen_problem.comm().get(),
1088  PetscObjectComposeFunction((PetscObject)mat,
1089  "formFunction",
1092 
1093  // It's also beneficial to be able to evaluate both A and B residuals at once
1094  LibmeshPetscCallA(eigen_problem.comm().get(),
1095  PetscObjectComposeFunction((PetscObject)mat,
1096  "formFunctionAB",
1098 
1099  // Users may choose to provide a custom measure of the norm of B (Bx for a linear system)
1100  if (eigen_problem.bxNormProvided())
1101  LibmeshPetscCallA(eigen_problem.comm().get(),
1102  PetscObjectComposeFunction((PetscObject)mat,
1103  "formNorm",
1105 
1106  // Finally we need to attach the "context" object, which is our EigenProblem, to the matrices so
1107  // that eventually when we get callbacks from SLEPc we can call methods on the EigenProblem
1108  PetscContainer container;
1109  LibmeshPetscCallA(eigen_problem.comm().get(),
1110  PetscContainerCreate(eigen_problem.comm().get(), &container));
1111  LibmeshPetscCallA(eigen_problem.comm().get(),
1112  PetscContainerSetPointer(container, &eigen_problem));
1113  LibmeshPetscCallA(
1114  eigen_problem.comm().get(),
1115  PetscObjectCompose((PetscObject)mat, "formJacobianCtx", (PetscObject)container));
1116  LibmeshPetscCallA(
1117  eigen_problem.comm().get(),
1118  PetscObjectCompose((PetscObject)mat, "formFunctionCtx", (PetscObject)container));
1119  if (eigen_problem.bxNormProvided())
1120  LibmeshPetscCallA(eigen_problem.comm().get(),
1121  PetscObjectCompose((PetscObject)mat, "formNormCtx", (PetscObject)container));
1122 
1123  LibmeshPetscCallA(eigen_problem.comm().get(), PetscContainerDestroy(&container));
1124 }
1125 
1126 PetscErrorCode
1127 mooseMatMult_Eigen(Mat mat, Vec x, Vec r)
1128 {
1130  void * ctx = nullptr;
1131  LibmeshPetscCallQ(MatShellGetContext(mat, &ctx));
1132 
1133  if (!ctx)
1134  mooseError("No context is set for shell matrix ");
1135 
1136  EigenProblem * eigen_problem = static_cast<EigenProblem *>(ctx);
1137  NonlinearEigenSystem & eigen_nl = eigen_problem->getCurrentNonlinearEigenSystem();
1138 
1139  evaluateResidual(*eigen_problem, x, r, eigen_nl.eigenVectorTag());
1140 
1141  if (eigen_problem->negativeSignEigenKernel())
1142  LibmeshPetscCallQ(VecScale(r, -1.));
1143 
1144  PetscFunctionReturn(PETSC_SUCCESS);
1145 }
1146 
1147 PetscErrorCode
1148 mooseMatMult_NonEigen(Mat mat, Vec x, Vec r)
1149 {
1151  void * ctx = nullptr;
1152  LibmeshPetscCallQ(MatShellGetContext(mat, &ctx));
1153 
1154  if (!ctx)
1155  mooseError("No context is set for shell matrix ");
1156 
1157  EigenProblem * eigen_problem = static_cast<EigenProblem *>(ctx);
1158  NonlinearEigenSystem & eigen_nl = eigen_problem->getCurrentNonlinearEigenSystem();
1159 
1160  evaluateResidual(*eigen_problem, x, r, eigen_nl.nonEigenVectorTag());
1161 
1162  PetscFunctionReturn(PETSC_SUCCESS);
1163 }
1164 
1165 void
1166 setOperationsForShellMat(EigenProblem & eigen_problem, Mat mat, bool eigen)
1167 {
1168  LibmeshPetscCallA(eigen_problem.comm().get(), MatShellSetContext(mat, &eigen_problem));
1169  LibmeshPetscCallA(eigen_problem.comm().get(),
1170  MatShellSetOperation(mat,
1171  MATOP_MULT,
1172  eigen ? (void (*)(void))mooseMatMult_Eigen
1173  : (void (*)(void))mooseMatMult_NonEigen));
1174 }
1175 
1176 PETSC_EXTERN PetscErrorCode
1178 {
1180 
1181  LibmeshPetscCallQ(PCRegister("moosepc", PCCreate_MoosePC));
1182 
1183  PetscFunctionReturn(PETSC_SUCCESS);
1184 }
1185 
1186 PETSC_EXTERN PetscErrorCode
1188 {
1190 
1191  pc->ops->view = PCView_MoosePC;
1192  pc->ops->destroy = PCDestroy_MoosePC;
1193  pc->ops->setup = PCSetUp_MoosePC;
1194  pc->ops->apply = PCApply_MoosePC;
1195 
1196  PetscFunctionReturn(PETSC_SUCCESS);
1197 }
1198 
1199 PetscErrorCode
1201 {
1203  /* We do not need to do anything right now, but later we may have some data we need to free here
1204  */
1205  PetscFunctionReturn(PETSC_SUCCESS);
1206 }
1207 
1208 PetscErrorCode
1209 PCView_MoosePC(PC /*pc*/, PetscViewer viewer)
1210 {
1211  PetscBool iascii;
1212 
1214  LibmeshPetscCallQ(PetscObjectTypeCompare((PetscObject)viewer, PETSCVIEWERASCII, &iascii));
1215  if (iascii)
1216  LibmeshPetscCallQ(PetscViewerASCIIPrintf(viewer, " %s\n", "moosepc"));
1217 
1218  PetscFunctionReturn(PETSC_SUCCESS);
1219 }
1220 
1221 PetscErrorCode
1222 PCApply_MoosePC(PC pc, Vec x, Vec y)
1223 {
1224  void * ctx;
1225  Mat Amat, Pmat;
1226  PetscContainer container;
1227 
1229  LibmeshPetscCallQ(PCGetOperators(pc, &Amat, &Pmat));
1231  PetscObjectQuery((PetscObject)Pmat, "formFunctionCtx", (PetscObject *)&container));
1232  if (container)
1233  LibmeshPetscCallQ(PetscContainerGetPointer(container, &ctx));
1234  else
1235  mooseError(" Can not find a context \n");
1236 
1237  EigenProblem * eigen_problem = static_cast<EigenProblem *>(ctx);
1238  NonlinearEigenSystem & nl_eigen = eigen_problem->getCurrentNonlinearEigenSystem();
1239  auto preconditioner = nl_eigen.preconditioner();
1240 
1241  if (!preconditioner)
1242  mooseError("There is no moose preconditioner in nonlinear eigen system \n");
1243 
1244  PetscVector<Number> x_vec(x, preconditioner->comm());
1245  PetscVector<Number> y_vec(y, preconditioner->comm());
1246 
1247  preconditioner->apply(x_vec, y_vec);
1248 
1249  PetscFunctionReturn(PETSC_SUCCESS);
1250 }
1251 
1252 PetscErrorCode
1254 {
1255  void * ctx;
1256  Mat Amat, Pmat;
1257  PetscContainer container;
1258 
1260  LibmeshPetscCallQ(PCGetOperators(pc, &Amat, &Pmat));
1262  PetscObjectQuery((PetscObject)Pmat, "formFunctionCtx", (PetscObject *)&container));
1263  if (container)
1264  LibmeshPetscCallQ(PetscContainerGetPointer(container, &ctx));
1265  else
1266  mooseError(" Can not find a context \n");
1267 
1268  EigenProblem * eigen_problem = static_cast<EigenProblem *>(ctx);
1269  NonlinearEigenSystem & nl_eigen = eigen_problem->getCurrentNonlinearEigenSystem();
1270  Preconditioner<Number> * preconditioner = nl_eigen.preconditioner();
1271 
1272  if (!preconditioner)
1273  mooseError("There is no moose preconditioner in nonlinear eigen system \n");
1274 
1275  if (!preconditioner->initialized())
1276  preconditioner->init();
1277 
1278  preconditioner->setup();
1279 
1280  PetscFunctionReturn(PETSC_SUCCESS);
1281 }
1282 
1283 PetscErrorCode
1285  PetscInt its,
1286  PetscInt max_it,
1287  PetscInt nconv,
1288  PetscInt nev,
1289  EPSConvergedReason * reason,
1290  void * ctx)
1291 {
1292  EigenProblem * eigen_problem = static_cast<EigenProblem *>(ctx);
1293 
1295  LibmeshPetscCallQ(EPSStoppingBasic(eps, its, max_it, nconv, nev, reason, NULL));
1296 
1297  // If we do free power iteration, we need to mark the solver as converged.
1298  // It is because SLEPc does not offer a way to copy unconverged solution.
1299  // If the solver is not marked as "converged", we have no way to get solution
1300  // from slepc. Note marking as "converged" has no side-effects at all for us.
1301  // If free power iteration is used as a stand-alone solver, we won't trigger
1302  // as "doFreePowerIteration()" is false.
1303  if (eigen_problem->doFreePowerIteration() && its == max_it && *reason <= 0)
1304  {
1305  *reason = EPS_CONVERGED_USER;
1306  eps->nconv = 1;
1307  }
1308  PetscFunctionReturn(PETSC_SUCCESS);
1309 }
1310 
1311 PetscErrorCode
1312 mooseSlepcEPSGetSNES(EPS eps, SNES * snes)
1313 {
1314  PetscBool same, nonlinear;
1315 
1317  LibmeshPetscCallQ(PetscObjectTypeCompare((PetscObject)eps, EPSPOWER, &same));
1318 
1319  if (!same)
1320  mooseError("It is not eps power, and there is no snes");
1321 
1322  LibmeshPetscCallQ(EPSPowerGetNonlinear(eps, &nonlinear));
1323 
1324  if (!nonlinear)
1325  mooseError("It is not a nonlinear eigen solver");
1326 
1327  LibmeshPetscCallQ(EPSPowerGetSNES(eps, snes));
1328 
1329  PetscFunctionReturn(PETSC_SUCCESS);
1330 }
1331 
1332 PetscErrorCode
1334 {
1335  SNES snes;
1336  const char * prefix = nullptr;
1337 
1340  // There is an extra "eps_power" in snes that users do not like it.
1341  // Let us remove that from snes.
1342  // Retrieve option prefix from EPS
1343  LibmeshPetscCallQ(PetscObjectGetOptionsPrefix((PetscObject)eps, &prefix));
1344  // Set option prefix to SNES
1345  LibmeshPetscCallQ(SNESSetOptionsPrefix(snes, prefix));
1346 
1347  PetscFunctionReturn(PETSC_SUCCESS);
1348 }
1349 
1350 PetscErrorCode
1352 {
1353  SNES snes;
1354  KSP ksp;
1355  PC pc;
1356 
1358  // Get SNES from EPS
1360  // Get KSP from SNES
1361  LibmeshPetscCallQ(SNESGetKSP(snes, &ksp));
1362  // Get PC from KSP
1363  LibmeshPetscCallQ(KSPGetPC(ksp, &pc));
1364  // Set PC type
1365  LibmeshPetscCallQ(PCSetType(pc, "moosepc"));
1366  PetscFunctionReturn(PETSC_SUCCESS);
1367 }
1368 
1369 PetscErrorCode
1371 {
1372  SNES snes;
1373  KSP ksp;
1374 
1376  // Get SNES from EPS
1378  // Get KSP from SNES
1379  LibmeshPetscCallQ(SNESGetKSP(snes, &ksp));
1380 
1382 
1384  PetscFunctionReturn(PETSC_SUCCESS);
1385 }
1386 
1387 PetscErrorCode
1389  PetscInt its,
1390  PetscInt /*nconv*/,
1391  PetscScalar * eigr,
1392  PetscScalar * eigi,
1393  PetscReal * /*errest*/,
1394  PetscInt /*nest*/,
1395  void * mctx)
1396 {
1397  ST st;
1398  PetscScalar eigenr, eigeni;
1399 
1401  EigenProblem * eigen_problem = static_cast<EigenProblem *>(mctx);
1402  auto & console = eigen_problem->console();
1403 
1404  auto inverse = eigen_problem->outputInverseEigenvalue();
1405  LibmeshPetscCallQ(EPSGetST(eps, &st));
1406  eigenr = eigr[0];
1407  eigeni = eigi[0];
1408  // Make the eigenvalue consistent with shift type
1409  LibmeshPetscCallQ(STBackTransform(st, 1, &eigenr, &eigeni));
1410 
1411  auto eigenvalue = inverse ? 1.0 / eigenr : eigenr;
1412 
1413  // The term "k-eigenvalue" is adopted from the neutronics community.
1414  console << " Iteration " << its << std::setprecision(10) << std::fixed
1415  << (inverse ? " k-eigenvalue = " : " eigenvalue = ") << eigenvalue << std::endl;
1416 
1417  PetscFunctionReturn(PETSC_SUCCESS);
1418 }
1419 
1420 } // namespace SlepcSupport
1421 } // namespace moose
1422 
1423 #endif // LIBMESH_HAVE_SLEPC
void setNewtonPetscOptions(SolverParams &solver_params, const InputParameters &params)
Definition: SlepcSupport.C:432
Nonlinear eigenvalue system to be solved.
bool constantMatrices() const
Whether or not matrices are constant.
Definition: EigenProblem.h:207
bool wereMatricesFormed() const
Whether or not constant matrices were already formed.
Definition: EigenProblem.h:217
Generalized Non-Hermitian.
Definition: MooseTypes.h:877
virtual void computeResidualTag(const NumericVector< Number > &soln, NumericVector< Number > &residual, TagID tag) override
Form a vector for all kernels and BCs with a given tag.
Definition: EigenProblem.C:295
unsigned int _solver_sys_num
Definition: SolverParams.h:36
Moose::PetscSupport::PetscOptions & getPetscOptions()
Retrieve a writable reference the PETSc options (used by PetscSupport)
const SparseMatrix< Number > & get_precond_matrix() const
int eps(unsigned int i, unsigned int j)
2D version
Generalized Hermitian indefinite.
Definition: MooseTypes.h:876
void addPetscOptionsFromCommandline()
Definition: PetscSupport.C:191
bool _eigen_matrix_free
Definition: SolverParams.h:27
Newton-based eigensolver with an assembled Jacobian matrix (fully coupled by default) ...
Definition: MooseTypes.h:862
Moose::EigenSolveType _eigen_solve_type
Definition: SolverParams.h:24
smallest magnitude
Definition: MooseTypes.h:888
const unsigned int invalid_uint
PetscErrorCode mooseSlepcEigenFormFunctionMFFD(void *ctx, Vec x, Vec r)
Function call for MFFD.
Definition: SlepcSupport.C:750
std::string _prefix
Definition: SolverParams.h:35
unsigned int TagID
Definition: MooseTypes.h:210
virtual std::size_t numNonlinearSystems() const override
bool outputInverseEigenvalue() const
Whether or not to output eigenvalue inverse.
Definition: EigenProblem.h:187
const SparseMatrix< Number > & get_matrix_A() const
void setDocString(const std::string &name, const std::string &doc)
Set the doc string of a parameter.
void moosePetscSNESFormMatricesTags(SNES, Vec x, std::vector< Mat > &eigen_mats, std::vector< SparseMatrix< Number > *> &all_dofs_mats, void *ctx, const std::set< TagID > &tags)
Definition: SlepcSupport.C:713
void mooseError(Args &&... args)
Emit an error message with the given stringified, concatenated args and terminate the application...
Definition: MooseError.h:333
PetscErrorCode mooseSlepcEPSSNESSetCustomizePC(EPS eps)
Attach a customized PC.
PetscErrorCode mooseSlepcStoppingTest(EPS eps, PetscInt its, PetscInt max_it, PetscInt nconv, PetscInt nev, EPSConvergedReason *reason, void *ctx)
A customized convergence checker.
std::vector< std::pair< R1, R2 > > get(const std::string &param1, const std::string &param2) const
Combine two vector parameters into a single vector of pairs.
TagID nonEigenVectorTag() const
Vector tag ID of left hand side.
void setCoupling(Moose::CouplingType type)
Set the coupling between variables TODO: allow user-defined coupling.
const SparseMatrix< Number > & get_matrix_B() const
void setEigenSolverOptions(SolverParams &solver_params, const InputParameters &params)
Definition: SlepcSupport.C:495
PetscErrorCode mooseMatMult_Eigen(Mat mat, Vec x, Vec y)
Implement MatMult via function evaluation for Bx.
The same as PJFNK except that matrix-vector multiplication is employed to replace residual evaluation...
Definition: MooseTypes.h:864
void mooseWarning(Args &&... args)
Emit a warning message with the given stringified, concatenated args.
Definition: MooseError.h:367
TagID eigenVectorTag() const
Vector tag ID of right hand side.
PetscErrorCode mooseSlepcEigenFormJacobianB(SNES snes, Vec x, Mat jac, Mat pc, void *ctx)
Form Jacobian matrix B.
Definition: SlepcSupport.C:862
void clearFreeNonlinearPowerIterations(const InputParameters &params)
Definition: SlepcSupport.C:420
T & set(const std::string &name, bool quiet_mode=false)
Returns a writable reference to the named parameters.
PetscFunctionBegin
virtual void execute(const ExecFlagType &exec_type) override
Convenience function for performing execution of MOOSE systems.
Definition: EigenProblem.C:165
PetscErrorCode mooseSlepcEigenFormFunctionAB(SNES snes, Vec x, Vec Ax, Vec Bx, void *ctx)
Form function residual Ax-Bx.
Definition: SlepcSupport.C:986
The main MOOSE class responsible for handling user-defined parameters in almost every MOOSE system...
bool negativeSignEigenKernel() const
A flag indicates if a negative sign is used in eigen kernels.
Definition: EigenProblem.h:57
const Parallel::Communicator & comm() const
void petscSetDefaultKSPNormType(FEProblemBase &problem, KSP ksp)
Set norm type.
Definition: PetscSupport.C:400
PetscErrorCode PCApply_MoosePC(PC pc, Vec x, Vec y)
Preconditioner application.
use whatever SLPEC has by default
Definition: MooseTypes.h:879
bool _customized_pc_for_eigen
Definition: SolverParams.h:29
target magnitude
Definition: MooseTypes.h:893
void inverse(const std::vector< std::vector< Real >> &m, std::vector< std::vector< Real >> &m_inv)
Inverse the dense square matrix m using LAPACK routines.
Definition: MatrixTools.C:23
PetscErrorCode PCDestroy_MoosePC(PC pc)
Destroy preconditioner.
PetscErrorCode mooseSlepcEigenFormFunctionA(SNES snes, Vec x, Vec r, void *ctx)
Form function residual Ax.
Definition: SlepcSupport.C:912
The following methods are specializations for using the libMesh::Parallel::packed_range_* routines fo...
void copy_super_to_sub(NumericVector< Number > &super, NumericVector< Number > &sub)
void setNEigenPairsRequired(unsigned int n_eigen_pairs)
Definition: EigenProblem.h:38
virtual bool initialized() const
Generalized Non-Hermitian with positive (semi-)definite B.
Definition: MooseTypes.h:878
void setEigenProblemOptions(SolverParams &solver_params, const MultiMooseEnum &dont_add_these_options)
Definition: SlepcSupport.C:299
Specialization of SubProblem for solving nonlinear equations plus auxiliary equations.
PetscErrorCode mooseMatMult_NonEigen(Mat mat, Vec x, Vec y)
Implement MatMult via function evaluation for Ax.
Power / Inverse / RQI.
Definition: MooseTypes.h:857
bool isNonlinearEigenvalueSolver(unsigned int eigen_sys_num) const
Definition: EigenProblem.C:676
InputParameters emptyInputParameters()
Krylov-Schur.
Definition: MooseTypes.h:859
void setSlepcEigenSolverTolerances(EigenProblem &eigen_problem, const SolverParams &solver_params, const InputParameters &params)
Control eigen solver tolerances via SLEPc options.
Definition: SlepcSupport.C:129
Real formNorm()
Form the Bx norm.
Definition: EigenProblem.C:693
virtual void swap(NumericVector< T > &v) override
InputParameters getSlepcValidParams(InputParameters &params)
Definition: SlepcSupport.C:38
virtual EPS getEPS()
Retrieve EPS (SLEPc eigen solver)
void libmesh_ignore(const Args &...)
bool _eigen_matrix_vector_mult
Definition: SolverParams.h:28
Preconditioned Jacobian-free Newton Krylov.
Definition: MooseTypes.h:863
PetscErrorCode PCSetUp_MoosePC(PC pc)
Setup preconditioner.
Moose::EigenProblemType _eigen_problem_type
Definition: SolverParams.h:25
virtual void zero()=0
void attachCallbacksToMat(EigenProblem &eigen_problem, Mat mat, bool eigen)
Attach call backs to mat.
void copy_sub_to_super(const NumericVector< Number > &sub, NumericVector< Number > &super)
void slepcSetOptions(EigenProblem &eigen_problem, SolverParams &solver_params, const InputParameters &params)
Push all SLEPc/PETSc options into SLEPc/PETSc side.
Definition: SlepcSupport.C:551
PETSC_EXTERN PetscErrorCode PCCreate_MoosePC(PC pc)
Create a preconditioner from moose side.
bool constJacobian() const
Returns _const_jacobian (whether a MOOSE object has specified that the Jacobian is the same as the pr...
PetscErrorCode mooseSlepcEigenFormNorm(SNES, Vec, PetscReal *norm, void *ctx)
virtual libMesh::EquationSystems & es() override
void computeResidualAB(const NumericVector< Number > &soln, NumericVector< Number > &residualA, NumericVector< Number > &residualB, TagID tagA, TagID tagB)
Form two vetors, where each is associated with one tag, through one element-loop. ...
Definition: EigenProblem.C:327
std::unique_ptr< NumericVector< Number > > solution
target imaginary
Definition: MooseTypes.h:895
This is a "smart" enum class intended to replace many of the shortcomings in the C++ enum type It sho...
Definition: MooseEnum.h:33
PetscErrorCode PCView_MoosePC(PC pc, PetscViewer viewer)
View preconditioner.
Moose::WhichEigenPairs _which_eigen_pairs
Definition: SolverParams.h:26
unsigned int number() const
Gets the number of this system.
Definition: SystemBase.C:1149
auto norm(const T &a) -> decltype(std::abs(a))
const ExecFlagType EXEC_LINEAR
Definition: Moose.C:31
void setNonlinearPowerOptions(SolverParams &solver_params)
Definition: SlepcSupport.C:475
std::string stringify(const T &t)
conversion to string
Definition: Conversion.h:64
void computeMatricesTags(const NumericVector< Number > &soln, const std::vector< SparseMatrix< Number > *> &jacobians, const std::set< TagID > &tags)
Form several matrices simultaneously.
Definition: EigenProblem.C:207
Nonlinear inverse power.
Definition: MooseTypes.h:861
void petscSetOptions(const PetscOptions &po, const SolverParams &solver_params, FEProblemBase *const problem=nullptr)
A function for setting the PETSc options in PETSc from the options supplied to MOOSE.
Definition: PetscSupport.C:230
NonlinearEigenSystem & getCurrentNonlinearEigenSystem()
Definition: EigenProblem.h:307
PetscErrorCode mooseSlepcEigenFormJacobianA(SNES snes, Vec x, Mat jac, Mat pc, void *ctx)
Form Jacobian matrix A.
Definition: SlepcSupport.C:776
PetscErrorCode mooseSlepcEigenFormFunctionB(SNES snes, Vec x, Vec r, void *ctx)
Form function residual Bx.
Definition: SlepcSupport.C:944
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.
Jacobi-Davidson.
Definition: MooseTypes.h:860
bool _precond_matrix_free
Definition: SolverParams.h:30
void moosePetscSNESFormFunction(SNES, Vec x, Vec r, void *ctx, TagID tag)
Definition: SlepcSupport.C:905
LibmeshPetscCallQ(DMShellGetContext(dm, &ctx))
const ConsoleStream & console() const
Return console handle.
Definition: Problem.h:48
void setSinglePetscOption(const std::string &name, const std::string &value="", FEProblemBase *const problem=nullptr)
A wrapper function for dealing with different versions of PetscOptionsSetValue.
DIE A HORRIBLE DEATH HERE typedef LIBMESH_DEFAULT_SCALAR_TYPE Real
PetscErrorCode mooseSlepcEPSGetSNES(EPS eps, SNES *snes)
Retrieve SNES from EPS.
bool onLinearSolver() const
Whether or not we are in a linear solver iteration.
Definition: EigenProblem.h:197
smallest imaginary
Definition: MooseTypes.h:892
void setFreeNonlinearPowerIterations(unsigned int free_power_iterations)
Set SLEPc/PETSc options to trigger free power iteration.
Definition: SlepcSupport.C:405
T & set(const std::string &)
libMesh::CondensedEigenSystem & sys()
Jacobian-free Newton Krylov.
Definition: MooseTypes.h:865
virtual SNES getSNES() override
Retrieve snes from slepc eigen solver.
TagID eigenMatrixTag() const
Matrix tag ID of right hand side.
largest imaginary
Definition: MooseTypes.h:891
IntRange< T > make_range(T beg, T end)
Non-Hermitian.
Definition: MooseTypes.h:874
PetscErrorCode mooseSlepcEPSSNESKSPSetPCSide(FEProblemBase &problem, EPS eps)
Allow users to specify PC side.
SolverParams & solverParams(unsigned int solver_sys_num=0)
Get the solver parameters.
std::unique_ptr< NumericVector< Number > > current_local_solution
void moosePetscSNESFormMatrixTag(SNES, Vec x, Mat eigen_mat, SparseMatrix< Number > &all_dofs_mat, void *ctx, TagID tag)
Definition: SlepcSupport.C:682
PetscErrorCode mooseSlepcEPSSNESSetUpOptionPrefix(EPS eps)
Get rid of prefix "-eps_power" for SNES, KSP, PC, etc.
all eigenvalues
Definition: MooseTypes.h:896
void addParam(const std::string &name, const S &value, const std::string &doc_string)
These methods add an optional parameter and a documentation string to the InputParameters object...
void setSinglePetscOptionIfAppropriate(const MultiMooseEnum &dont_add_these_options, const std::string &name, const std::string &value="", FEProblemBase *const problem=nullptr)
Same as setSinglePetscOption, but does not set the option if it doesn&#39;t make sense for the current si...
void * ctx
void setOperationsForShellMat(EigenProblem &eigen_problem, Mat mat, bool eigen)
Set operations to shell mat.
void setEigenProblemSolverParams(EigenProblem &eigen_problem, const InputParameters &params)
Retrieve eigen problem params from &#39;params&#39;, and then set these params into SolverParams.
Definition: SlepcSupport.C:213
TagID nonEigenMatrixTag() const
Matrix tag ID of left hand side.
MOOSE now contains C++17 code, so give a reasonable error message stating what the user can do to add...
This is a "smart" enum class intended to replace many of the shortcomings in the C++ enum type...
use whatever we have in SLEPC
Definition: MooseTypes.h:897
const int subspace_factor
Definition: SlepcSupport.C:35
PetscErrorCode mooseEPSFormMatrices(EigenProblem &eigen_problem, EPS eps, Vec x, void *ctx)
Definition: SlepcSupport.C:574
Problem for solving eigenvalue problems.
Definition: EigenProblem.h:21
void storeSolveType(FEProblemBase &fe_problem, const InputParameters &params)
Set solve type into eigen problem (solverParams)
Definition: SlepcSupport.C:287
largest magnitude
Definition: MooseTypes.h:887
MultiMooseEnum dont_add_these_options
Flags to explicitly not set, even if they are specified programmatically.
Definition: PetscSupport.h:59
PETSC_EXTERN PetscErrorCode registerPCToPETSc()
Let PETSc know there is a preconditioner.
bool bxNormProvided() const
Whether a Bx norm postprocessor has been provided.
Definition: EigenProblem.h:232
void setWhichEigenPairsOptions(SolverParams &solver_params, const MultiMooseEnum &dont_add_these_options)
Definition: SlepcSupport.C:342
PetscFunctionReturn(LIBMESH_PETSC_SUCCESS)
virtual void setup()
const DofMap & get_dof_map() const
bool doFreePowerIteration() const
Whether or not we are doing free power iteration.
Definition: EigenProblem.h:79
libMesh::Preconditioner< Number > * preconditioner() const
auto index_range(const T &sizable)
Generalized Hermitian.
Definition: MooseTypes.h:875
void petscSetDefaultPCSide(FEProblemBase &problem, KSP ksp)
Setup which side we want to apply preconditioner.
Definition: PetscSupport.C:411
virtual void computeJacobianTag(const NumericVector< Number > &soln, SparseMatrix< Number > &jacobian, TagID tag) override
Form a Jacobian matrix for all kernels and BCs with a given tag.
Definition: EigenProblem.C:176
InputParameters getSlepcEigenProblemValidParams()
Retrieve valid params that allow users to specify eigen problem configuration.
Definition: SlepcSupport.C:67
void addParamNamesToGroup(const std::string &space_delim_names, const std::string group_name)
This method takes a space delimited list of parameter names and adds them to the specified group name...
bool isParamValid(const std::string &name) const
This method returns parameters that have been initialized in one fashion or another, i.e.