Line data Source code
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 :
37 : InputParameters
38 15437 : getSlepcValidParams(InputParameters & params)
39 : {
40 : MooseEnum solve_type("POWER ARNOLDI KRYLOVSCHUR JACOBI_DAVIDSON "
41 : "NONLINEAR_POWER NEWTON PJFNK PJFNKMO JFNK",
42 61748 : "PJFNK");
43 30874 : params.set<MooseEnum>("solve_type") = solve_type;
44 :
45 61748 : 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 15437 : params.set<Real>("l_tol") = 1e-2;
62 :
63 30874 : return params;
64 15437 : }
65 :
66 : InputParameters
67 15437 : getSlepcEigenProblemValidParams()
68 : {
69 15437 : InputParameters params = emptyInputParameters();
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 61748 : "GEN_NON_HERMITIAN");
75 61748 : 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 61748 : "TARGET_IMAGINARY ALL_EIGENVALUES SLEPC_DEFAULT");
91 61748 : 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 61748 : params.addParam<unsigned int>("n_eigen_pairs", 1, "The number of eigen pairs");
107 61748 : params.addParam<unsigned int>("n_basis_vectors", 3, "The dimension of eigen subspaces");
108 :
109 61748 : params.addParam<Real>("eigen_tol", 1.0e-4, "Relative Tolerance for Eigen Solver");
110 61748 : params.addParam<unsigned int>("eigen_max_its", 10000, "Max Iterations for Eigen Solver");
111 :
112 61748 : params.addParam<Real>("l_abs_tol", 1e-50, "Absolute Tolerances for Linear Solver");
113 :
114 61748 : params.addParam<unsigned int>("free_power_iterations", 4, "The number of free power iterations");
115 :
116 46311 : params.addParam<unsigned int>(
117 30874 : "extra_power_iterations", 0, "The number of extra free power iterations");
118 :
119 61748 : 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 46311 : params.addParamNamesToGroup("l_abs_tol", "Linear solver");
124 :
125 30874 : return params;
126 15437 : }
127 :
128 : void
129 558 : setSlepcEigenSolverTolerances(EigenProblem & eigen_problem,
130 : const SolverParams & solver_params,
131 : const InputParameters & params)
132 : {
133 558 : const auto & dont_add_these_options = eigen_problem.getPetscOptions().dont_add_these_options;
134 558 : 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 :
139 558 : Moose::PetscSupport::setSinglePetscOptionIfAppropriate(dont_add_these_options,
140 1116 : prefix_with_dash + "eps_tol",
141 1116 : stringify(params.get<Real>("eigen_tol")));
142 :
143 558 : Moose::PetscSupport::setSinglePetscOptionIfAppropriate(
144 : dont_add_these_options,
145 1116 : prefix_with_dash + "eps_max_it",
146 1116 : 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 558 : if (eigen_problem.isNonlinearEigenvalueSolver(solver_params._solver_sys_num))
151 : {
152 : // nonlinear solver tolerances
153 493 : Moose::PetscSupport::setSinglePetscOptionIfAppropriate(
154 : dont_add_these_options,
155 986 : prefix_with_dash + "snes_max_it",
156 986 : stringify(params.get<unsigned int>("nl_max_its")));
157 :
158 493 : Moose::PetscSupport::setSinglePetscOptionIfAppropriate(
159 : dont_add_these_options,
160 986 : prefix_with_dash + "snes_max_funcs",
161 986 : stringify(params.get<unsigned int>("nl_max_funcs")));
162 :
163 493 : Moose::PetscSupport::setSinglePetscOptionIfAppropriate(
164 : dont_add_these_options,
165 986 : prefix_with_dash + "snes_atol",
166 986 : stringify(params.get<Real>("nl_abs_tol")));
167 :
168 493 : Moose::PetscSupport::setSinglePetscOptionIfAppropriate(
169 : dont_add_these_options,
170 986 : prefix_with_dash + "snes_rtol",
171 986 : stringify(params.get<Real>("nl_rel_tol")));
172 :
173 493 : Moose::PetscSupport::setSinglePetscOptionIfAppropriate(
174 : dont_add_these_options,
175 986 : prefix_with_dash + "snes_stol",
176 986 : stringify(params.get<Real>("nl_rel_step_tol")));
177 :
178 : // linear solver
179 493 : Moose::PetscSupport::setSinglePetscOptionIfAppropriate(
180 : dont_add_these_options,
181 986 : prefix_with_dash + "ksp_max_it",
182 986 : stringify(params.get<unsigned int>("l_max_its")));
183 :
184 493 : Moose::PetscSupport::setSinglePetscOptionIfAppropriate(dont_add_these_options,
185 986 : prefix_with_dash + "ksp_rtol",
186 986 : stringify(params.get<Real>("l_tol")));
187 :
188 493 : Moose::PetscSupport::setSinglePetscOptionIfAppropriate(
189 : dont_add_these_options,
190 986 : prefix_with_dash + "ksp_atol",
191 986 : stringify(params.get<Real>("l_abs_tol")));
192 : }
193 : else
194 : { // linear eigenvalue problem
195 : // linear solver
196 65 : Moose::PetscSupport::setSinglePetscOptionIfAppropriate(
197 : dont_add_these_options,
198 130 : prefix_with_dash + "st_ksp_max_it",
199 130 : stringify(params.get<unsigned int>("l_max_its")));
200 :
201 65 : Moose::PetscSupport::setSinglePetscOptionIfAppropriate(dont_add_these_options,
202 130 : prefix_with_dash + "st_ksp_rtol",
203 130 : stringify(params.get<Real>("l_tol")));
204 :
205 65 : Moose::PetscSupport::setSinglePetscOptionIfAppropriate(
206 : dont_add_these_options,
207 130 : prefix_with_dash + "st_ksp_atol",
208 130 : stringify(params.get<Real>("l_abs_tol")));
209 : }
210 558 : }
211 :
212 : void
213 588 : setEigenProblemSolverParams(EigenProblem & eigen_problem, const InputParameters & params)
214 : {
215 1176 : for (const auto i : make_range(eigen_problem.numNonlinearSystems()))
216 : {
217 588 : const std::string & eigen_problem_type = params.get<MooseEnum>("eigen_problem_type");
218 588 : if (!eigen_problem_type.empty())
219 588 : eigen_problem.solverParams(i)._eigen_problem_type =
220 588 : Moose::stringToEnum<Moose::EigenProblemType>(eigen_problem_type);
221 : else
222 0 : mooseError("Have to specify a valid eigen problem type");
223 :
224 588 : const std::string & which_eigen_pairs = params.get<MooseEnum>("which_eigen_pairs");
225 588 : if (!which_eigen_pairs.empty())
226 77 : eigen_problem.solverParams(i)._which_eigen_pairs =
227 77 : 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 588 : unsigned int n_eigen_pairs = params.get<unsigned int>("n_eigen_pairs");
234 588 : unsigned int n_basis_vectors = params.get<unsigned int>("n_basis_vectors");
235 :
236 588 : eigen_problem.setNEigenPairsRequired(n_eigen_pairs);
237 :
238 1176 : 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 588 : if (subspace_factor * n_eigen_pairs > n_basis_vectors)
242 : {
243 0 : n_basis_vectors = subspace_factor * n_eigen_pairs;
244 0 : mooseWarning(
245 : "Number of subspaces in Eigensolver is changed by moose because the value you set "
246 : "is too small");
247 : }
248 :
249 1176 : eigen_problem.es().parameters.set<unsigned int>("basis vectors") = n_basis_vectors;
250 :
251 : // Operators A and B are formed as shell matrices
252 588 : eigen_problem.solverParams(i)._eigen_matrix_free = params.get<bool>("matrix_free");
253 :
254 : // Preconditioning is formed as a shell matrix
255 588 : eigen_problem.solverParams(i)._precond_matrix_free = params.get<bool>("precond_matrix_free");
256 :
257 588 : if (params.get<MooseEnum>("solve_type") == "PJFNK")
258 : {
259 320 : eigen_problem.solverParams(i)._eigen_matrix_free = true;
260 : }
261 588 : if (params.get<MooseEnum>("solve_type") == "JFNK")
262 : {
263 40 : eigen_problem.solverParams(i)._eigen_matrix_free = true;
264 40 : eigen_problem.solverParams(i)._precond_matrix_free = true;
265 : }
266 : // We need matrices so that we can implement residual evaluations
267 588 : if (params.get<MooseEnum>("solve_type") == "PJFNKMO")
268 : {
269 102 : eigen_problem.solverParams(i)._eigen_matrix_free = true;
270 102 : eigen_problem.solverParams(i)._precond_matrix_free = false;
271 102 : 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 102 : eigen_problem.setCoupling(Moose::COUPLING_FULL);
275 : }
276 588 : }
277 :
278 588 : eigen_problem.constantMatrices(params.get<bool>("constant_matrices"));
279 :
280 588 : if (eigen_problem.constantMatrices() && params.get<MooseEnum>("solve_type") != "PJFNKMO")
281 : {
282 4 : mooseError("constant_matrices flag is only valid for solve type: PJFNKMO");
283 : }
284 584 : }
285 :
286 : void
287 588 : storeSolveType(FEProblemBase & fe_problem, const InputParameters & params)
288 : {
289 588 : if (!(dynamic_cast<EigenProblem *>(&fe_problem)))
290 0 : return;
291 :
292 1176 : if (params.isParamValid("solve_type"))
293 1176 : for (const auto i : make_range(fe_problem.numNonlinearSystems()))
294 1176 : fe_problem.solverParams(i)._eigen_solve_type =
295 588 : Moose::stringToEnum<Moose::EigenSolveType>(params.get<MooseEnum>("solve_type"));
296 : }
297 :
298 : void
299 558 : setEigenProblemOptions(SolverParams & solver_params, const MultiMooseEnum & dont_add_these_options)
300 : {
301 558 : switch (solver_params._eigen_problem_type)
302 : {
303 26 : case Moose::EPT_HERMITIAN:
304 78 : Moose::PetscSupport::setSinglePetscOptionIfAppropriate(dont_add_these_options,
305 : "-eps_hermitian");
306 26 : break;
307 :
308 13 : case Moose::EPT_NON_HERMITIAN:
309 39 : Moose::PetscSupport::setSinglePetscOptionIfAppropriate(dont_add_these_options,
310 : "-eps_non_hermitian");
311 13 : break;
312 :
313 0 : case Moose::EPT_GEN_HERMITIAN:
314 0 : Moose::PetscSupport::setSinglePetscOptionIfAppropriate(dont_add_these_options,
315 : "-eps_gen_hermitian");
316 0 : break;
317 :
318 0 : case Moose::EPT_GEN_INDEFINITE:
319 0 : Moose::PetscSupport::setSinglePetscOptionIfAppropriate(dont_add_these_options,
320 : "-eps_gen_indefinite");
321 0 : break;
322 :
323 519 : case Moose::EPT_GEN_NON_HERMITIAN:
324 1557 : Moose::PetscSupport::setSinglePetscOptionIfAppropriate(dont_add_these_options,
325 : "-eps_gen_non_hermitian");
326 519 : break;
327 :
328 0 : case Moose::EPT_POS_GEN_NON_HERMITIAN:
329 0 : Moose::PetscSupport::setSinglePetscOptionIfAppropriate(dont_add_these_options,
330 : "-eps_pos_gen_non_hermitian");
331 0 : break;
332 :
333 0 : case Moose::EPT_SLEPC_DEFAULT:
334 0 : break;
335 :
336 0 : default:
337 0 : mooseError("Unknown eigen solver type \n");
338 : }
339 558 : }
340 :
341 : void
342 558 : setWhichEigenPairsOptions(SolverParams & solver_params,
343 : const MultiMooseEnum & dont_add_these_options)
344 : {
345 558 : switch (solver_params._which_eigen_pairs)
346 : {
347 39 : case Moose::WEP_LARGEST_MAGNITUDE:
348 117 : Moose::PetscSupport::setSinglePetscOptionIfAppropriate(dont_add_these_options,
349 : "-eps_largest_magnitude");
350 39 : break;
351 :
352 26 : case Moose::WEP_SMALLEST_MAGNITUDE:
353 78 : Moose::PetscSupport::setSinglePetscOptionIfAppropriate(dont_add_these_options,
354 : "-eps_smallest_magnitude");
355 26 : break;
356 :
357 0 : case Moose::WEP_LARGEST_REAL:
358 0 : Moose::PetscSupport::setSinglePetscOptionIfAppropriate(dont_add_these_options,
359 : "-eps_largest_real");
360 0 : break;
361 :
362 0 : case Moose::WEP_SMALLEST_REAL:
363 0 : Moose::PetscSupport::setSinglePetscOptionIfAppropriate(dont_add_these_options,
364 : "-eps_smallest_real");
365 0 : break;
366 :
367 0 : case Moose::WEP_LARGEST_IMAGINARY:
368 0 : Moose::PetscSupport::setSinglePetscOptionIfAppropriate(dont_add_these_options,
369 : "-eps_largest_imaginary");
370 0 : break;
371 :
372 0 : case Moose::WEP_SMALLEST_IMAGINARY:
373 0 : Moose::PetscSupport::setSinglePetscOptionIfAppropriate(dont_add_these_options,
374 : "-eps_smallest_imaginary");
375 0 : break;
376 :
377 0 : case Moose::WEP_TARGET_MAGNITUDE:
378 0 : Moose::PetscSupport::setSinglePetscOptionIfAppropriate(dont_add_these_options,
379 : "-eps_target_magnitude");
380 0 : break;
381 :
382 0 : case Moose::WEP_TARGET_REAL:
383 0 : Moose::PetscSupport::setSinglePetscOptionIfAppropriate(dont_add_these_options,
384 : "-eps_target_real");
385 0 : break;
386 :
387 0 : case Moose::WEP_TARGET_IMAGINARY:
388 0 : Moose::PetscSupport::setSinglePetscOptionIfAppropriate(dont_add_these_options,
389 : "-eps_target_imaginary");
390 0 : break;
391 :
392 0 : case Moose::WEP_ALL_EIGENVALUES:
393 0 : Moose::PetscSupport::setSinglePetscOptionIfAppropriate(dont_add_these_options, "-eps_all");
394 0 : break;
395 :
396 493 : case Moose::WEP_SLEPC_DEFAULT:
397 493 : break;
398 :
399 0 : default:
400 0 : mooseError("Unknown type of WhichEigenPairs \n");
401 : }
402 558 : }
403 :
404 : void
405 558 : setFreeNonlinearPowerIterations(unsigned int free_power_iterations)
406 : {
407 2232 : Moose::PetscSupport::setSinglePetscOption("-eps_power_update", "0");
408 2232 : 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 1674 : Moose::PetscSupport::setSinglePetscOption("-snes_rtol", "0.99999999999");
413 1674 : 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 1674 : Moose::PetscSupport::setSinglePetscOption("-eps_tol", "1e-50");
417 558 : }
418 :
419 : void
420 558 : clearFreeNonlinearPowerIterations(const InputParameters & params)
421 : {
422 2232 : Moose::PetscSupport::setSinglePetscOption("-eps_power_update", "1");
423 1674 : Moose::PetscSupport::setSinglePetscOption("-eps_max_it", "1");
424 558 : Moose::PetscSupport::setSinglePetscOption("-snes_max_it",
425 1116 : stringify(params.get<unsigned int>("nl_max_its")));
426 558 : Moose::PetscSupport::setSinglePetscOption("-snes_rtol",
427 1116 : stringify(params.get<Real>("nl_rel_tol")));
428 1674 : Moose::PetscSupport::setSinglePetscOption("-eps_tol", stringify(params.get<Real>("eigen_tol")));
429 558 : }
430 :
431 : void
432 480 : 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 480 : bool initial_power = params.get<bool>("_newton_inverse_power");
437 :
438 1920 : Moose::PetscSupport::setSinglePetscOption("-eps_type", "power");
439 1920 : Moose::PetscSupport::setSinglePetscOption("-eps_power_nonlinear", "1");
440 1920 : 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 1440 : Moose::PetscSupport::setSinglePetscOption("-eps_max_it", "1");
444 480 : if (initial_power)
445 : {
446 0 : Moose::PetscSupport::setSinglePetscOption("-init_eps_power_snes_max_it", "1");
447 0 : Moose::PetscSupport::setSinglePetscOption("-init_eps_power_ksp_rtol", "1e-2");
448 0 : Moose::PetscSupport::setSinglePetscOption(
449 0 : "-init_eps_max_it", stringify(params.get<unsigned int>("free_power_iterations")));
450 : }
451 1440 : Moose::PetscSupport::setSinglePetscOption("-eps_target_magnitude", "");
452 480 : if (solver_params._eigen_matrix_free)
453 : {
454 1401 : Moose::PetscSupport::setSinglePetscOption("-snes_mf_operator", "1");
455 467 : if (initial_power)
456 0 : Moose::PetscSupport::setSinglePetscOption("-init_eps_power_snes_mf_operator", "1");
457 : }
458 : else
459 : {
460 39 : Moose::PetscSupport::setSinglePetscOption("-snes_mf_operator", "0");
461 13 : if (initial_power)
462 0 : 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 480 : }
473 :
474 : void
475 13 : setNonlinearPowerOptions(SolverParams & solver_params)
476 : {
477 : #if !SLEPC_VERSION_LESS_THAN(3, 8, 0) || !PETSC_VERSION_RELEASE
478 52 : Moose::PetscSupport::setSinglePetscOption("-eps_type", "power");
479 52 : Moose::PetscSupport::setSinglePetscOption("-eps_power_nonlinear", "1");
480 39 : Moose::PetscSupport::setSinglePetscOption("-eps_target_magnitude", "");
481 13 : if (solver_params._eigen_matrix_free)
482 52 : Moose::PetscSupport::setSinglePetscOption("-snes_mf_operator", "1");
483 : else
484 0 : 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 13 : }
493 :
494 : void
495 558 : setEigenSolverOptions(SolverParams & solver_params, const InputParameters & params)
496 : {
497 : // Avoid unused variable warnings when you have SLEPc but not PETSc-dev.
498 558 : libmesh_ignore(params);
499 :
500 558 : switch (solver_params._eigen_solve_type)
501 : {
502 0 : case Moose::EST_POWER:
503 0 : Moose::PetscSupport::setSinglePetscOption("-eps_type", "power");
504 0 : break;
505 :
506 0 : case Moose::EST_ARNOLDI:
507 0 : Moose::PetscSupport::setSinglePetscOption("-eps_type", "arnoldi");
508 0 : break;
509 :
510 39 : case Moose::EST_KRYLOVSCHUR:
511 117 : Moose::PetscSupport::setSinglePetscOption("-eps_type", "krylovschur");
512 39 : break;
513 :
514 26 : case Moose::EST_JACOBI_DAVIDSON:
515 78 : Moose::PetscSupport::setSinglePetscOption("-eps_type", "jd");
516 26 : break;
517 :
518 13 : case Moose::EST_NONLINEAR_POWER:
519 13 : setNonlinearPowerOptions(solver_params);
520 13 : break;
521 :
522 36 : case Moose::EST_NEWTON:
523 36 : setNewtonPetscOptions(solver_params, params);
524 36 : break;
525 :
526 306 : case Moose::EST_PJFNK:
527 306 : solver_params._eigen_matrix_free = true;
528 306 : solver_params._customized_pc_for_eigen = false;
529 306 : setNewtonPetscOptions(solver_params, params);
530 306 : break;
531 :
532 40 : case Moose::EST_JFNK:
533 40 : solver_params._eigen_matrix_free = true;
534 40 : solver_params._customized_pc_for_eigen = true;
535 40 : setNewtonPetscOptions(solver_params, params);
536 40 : break;
537 :
538 98 : case Moose::EST_PJFNKMO:
539 98 : solver_params._eigen_matrix_free = true;
540 98 : solver_params._customized_pc_for_eigen = false;
541 98 : solver_params._eigen_matrix_vector_mult = true;
542 98 : setNewtonPetscOptions(solver_params, params);
543 98 : break;
544 :
545 0 : default:
546 0 : mooseError("Unknown eigen solver type \n");
547 : }
548 558 : }
549 :
550 : void
551 558 : slepcSetOptions(EigenProblem & eigen_problem,
552 : SolverParams & solver_params,
553 : const InputParameters & params)
554 : {
555 558 : const auto & dont_add_these_options = eigen_problem.getPetscOptions().dont_add_these_options;
556 :
557 558 : Moose::PetscSupport::petscSetOptions(
558 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 558 : setSlepcEigenSolverTolerances(eigen_problem, solver_params, params);
562 558 : setEigenSolverOptions(solver_params, params);
563 : // when Bx norm postprocessor is provided, we switch off the sign normalization
564 558 : if (eigen_problem.bxNormProvided())
565 165 : Moose::PetscSupport::setSinglePetscOptionIfAppropriate(
566 : dont_add_these_options, "-eps_power_sign_normalization", "0", &eigen_problem);
567 558 : setEigenProblemOptions(solver_params, eigen_problem.getPetscOptions().dont_add_these_options);
568 558 : setWhichEigenPairsOptions(solver_params, eigen_problem.getPetscOptions().dont_add_these_options);
569 558 : Moose::PetscSupport::addPetscOptionsFromCommandline();
570 558 : }
571 :
572 : // For matrices A and B
573 : PetscErrorCode
574 6160 : mooseEPSFormMatrices(EigenProblem & eigen_problem, EPS eps, Vec x, void * ctx)
575 : {
576 : ST st;
577 : Mat A, B;
578 : PetscBool aisshell, bisshell;
579 : PetscFunctionBegin;
580 :
581 6160 : if (eigen_problem.constantMatrices() && eigen_problem.wereMatricesFormed())
582 3368 : PetscFunctionReturn(PETSC_SUCCESS);
583 :
584 2792 : 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 2412 : PetscFunctionReturn(PETSC_SUCCESS);
589 :
590 380 : NonlinearEigenSystem & eigen_nl = eigen_problem.getCurrentNonlinearEigenSystem();
591 380 : auto & sys = eigen_nl.sys();
592 380 : SNES snes = eigen_nl.getSNES();
593 : // Rest ST state so that we can retrieve matrices
594 380 : LibmeshPetscCallQ(EPSGetST(eps, &st));
595 380 : LibmeshPetscCallQ(STResetMatrixState(st));
596 380 : LibmeshPetscCallQ(EPSGetOperators(eps, &A, &B));
597 380 : LibmeshPetscCallQ(PetscObjectTypeCompare((PetscObject)A, MATSHELL, &aisshell));
598 380 : LibmeshPetscCallQ(PetscObjectTypeCompare((PetscObject)B, MATSHELL, &bisshell));
599 380 : if (aisshell || bisshell)
600 : {
601 0 : 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 760 : std::vector<Mat> mats = {A, B};
607 760 : std::vector<SparseMatrix<Number> *> libmesh_mats = {&sys.get_matrix_A(), &sys.get_matrix_B()};
608 760 : moosePetscSNESFormMatricesTags(
609 380 : snes, x, mats, libmesh_mats, ctx, {eigen_nl.nonEigenMatrixTag(), eigen_nl.eigenMatrixTag()});
610 380 : eigen_problem.wereMatricesFormed(true);
611 380 : PetscFunctionReturn(PETSC_SUCCESS);
612 380 : }
613 :
614 : namespace
615 : {
616 : void
617 35238 : updateCurrentLocalSolution(CondensedEigenSystem & sys, Vec x)
618 : {
619 35238 : auto & dof_map = sys.get_dof_map();
620 :
621 35238 : PetscVector<Number> X_global(x, sys.comm());
622 :
623 35238 : if (dof_map.n_constrained_dofs())
624 : {
625 984 : sys.copy_sub_to_super(X_global, *sys.solution);
626 : // Set the constrained dof values
627 984 : dof_map.enforce_constraints_exactly(sys);
628 984 : sys.update();
629 : }
630 : else
631 : {
632 34254 : 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 34254 : X_global.swap(X_sys);
639 34254 : sys.update();
640 34254 : X_global.swap(X_sys);
641 : }
642 35238 : }
643 :
644 : std::unique_ptr<NumericVector<Number>>
645 44194 : createWrappedResidual(CondensedEigenSystem & sys, Vec r)
646 : {
647 44194 : auto & dof_map = sys.get_dof_map();
648 :
649 44194 : if (dof_map.n_constrained_dofs())
650 1296 : return sys.solution->zero_clone();
651 : else
652 : {
653 42898 : auto R = std::make_unique<PetscVector<Number>>(r, sys.comm());
654 42898 : R->zero();
655 42898 : return R;
656 42898 : }
657 : }
658 :
659 : void
660 17694 : evaluateResidual(EigenProblem & eigen_problem, Vec x, Vec r, TagID tag)
661 : {
662 17694 : auto & nl = eigen_problem.getCurrentNonlinearEigenSystem();
663 17694 : auto & sys = nl.sys();
664 17694 : auto & dof_map = sys.get_dof_map();
665 :
666 17694 : updateCurrentLocalSolution(sys, x);
667 17694 : auto R = createWrappedResidual(sys, r);
668 :
669 17694 : eigen_problem.computeResidualTag(*sys.current_local_solution.get(), *R, tag);
670 :
671 17694 : R->close();
672 :
673 17694 : if (dof_map.n_constrained_dofs())
674 : {
675 456 : PetscVector<Number> sub_r(r, sys.comm());
676 456 : sys.copy_super_to_sub(*R, sub_r);
677 456 : }
678 17694 : }
679 : }
680 :
681 : void
682 3650 : moosePetscSNESFormMatrixTag(
683 : SNES /*snes*/, Vec x, Mat eigen_mat, SparseMatrix<Number> & all_dofs_mat, void * ctx, TagID tag)
684 : {
685 3650 : EigenProblem * eigen_problem = static_cast<EigenProblem *>(ctx);
686 3650 : auto & nl = eigen_problem->getCurrentNonlinearEigenSystem();
687 3650 : auto & sys = nl.sys();
688 3650 : 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 3650 : updateCurrentLocalSolution(sys, x);
699 :
700 3650 : if (!eigen_problem->constJacobian())
701 3650 : all_dofs_mat.zero();
702 :
703 3650 : eigen_problem->computeJacobianTag(*sys.current_local_solution.get(), all_dofs_mat, tag);
704 :
705 3650 : if (dof_map.n_constrained_dofs())
706 : {
707 96 : PetscMatrix<Number> wrapped_eigen_mat(eigen_mat, sys.comm());
708 96 : sys.copy_super_to_sub(all_dofs_mat, wrapped_eigen_mat);
709 96 : }
710 3650 : }
711 :
712 : void
713 380 : moosePetscSNESFormMatricesTags(SNES /*snes*/,
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 380 : EigenProblem * eigen_problem = static_cast<EigenProblem *>(ctx);
721 380 : auto & nl = eigen_problem->getCurrentNonlinearEigenSystem();
722 380 : auto & sys = nl.sys();
723 380 : 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 380 : updateCurrentLocalSolution(sys, x);
734 :
735 1140 : for (auto * const all_dofs_mat : all_dofs_mats)
736 760 : if (!eigen_problem->constJacobian())
737 760 : all_dofs_mat->zero();
738 :
739 380 : eigen_problem->computeMatricesTags(*sys.current_local_solution.get(), all_dofs_mats, tags);
740 :
741 380 : if (dof_map.n_constrained_dofs())
742 36 : for (const auto i : index_range(eigen_mats))
743 : {
744 24 : PetscMatrix<Number> wrapped_eigen_mat(eigen_mats[i], sys.comm());
745 24 : sys.copy_super_to_sub(*all_dofs_mats[i], wrapped_eigen_mat);
746 24 : }
747 380 : }
748 :
749 : PetscErrorCode
750 4636 : mooseSlepcEigenFormFunctionMFFD(void * ctx, Vec x, Vec r)
751 : {
752 : PetscErrorCode (*func)(SNES, Vec, Vec, void *);
753 : void * fctx;
754 : PetscFunctionBegin;
755 :
756 4636 : EigenProblem * eigen_problem = static_cast<EigenProblem *>(ctx);
757 4636 : NonlinearEigenSystem & eigen_nl = eigen_problem->getCurrentNonlinearEigenSystem();
758 4636 : SNES snes = eigen_nl.getSNES();
759 :
760 4636 : eigen_problem->onLinearSolver(true);
761 :
762 4636 : LibmeshPetscCallQ(SNESGetFunction(snes, NULL, &func, &fctx));
763 4636 : if (fctx != ctx)
764 : {
765 0 : SETERRQ(
766 : PetscObjectComm((PetscObject)snes), PETSC_ERR_ARG_INCOMP, "Contexts are not consistent \n");
767 : }
768 4636 : LibmeshPetscCallQ((*func)(snes, x, r, ctx));
769 :
770 4636 : eigen_problem->onLinearSolver(false);
771 :
772 4636 : PetscFunctionReturn(PETSC_SUCCESS);
773 : }
774 :
775 : PetscErrorCode
776 4628 : mooseSlepcEigenFormJacobianA(SNES snes, Vec x, Mat jac, Mat pc, void * ctx)
777 : {
778 : PetscBool jisshell, pisshell;
779 : PetscBool jismffd;
780 :
781 : PetscFunctionBegin;
782 :
783 4628 : EigenProblem * eigen_problem = static_cast<EigenProblem *>(ctx);
784 4628 : NonlinearEigenSystem & eigen_nl = eigen_problem->getCurrentNonlinearEigenSystem();
785 4628 : auto & sys = eigen_nl.sys();
786 :
787 : // If both jacobian and preconditioning are shell matrices,
788 : // and then assemble them and return
789 4628 : LibmeshPetscCallQ(PetscObjectTypeCompare((PetscObject)jac, MATSHELL, &jisshell));
790 4628 : LibmeshPetscCallQ(PetscObjectTypeCompare((PetscObject)jac, MATMFFD, &jismffd));
791 :
792 4628 : if (jismffd && eigen_problem->solverParams(eigen_nl.number())._eigen_matrix_vector_mult)
793 : {
794 688 : LibmeshPetscCallQ(
795 : MatMFFDSetFunction(jac, Moose::SlepcSupport::mooseSlepcEigenFormFunctionMFFD, ctx));
796 :
797 688 : EPS eps = eigen_nl.getEPS();
798 :
799 688 : LibmeshPetscCallQ(mooseEPSFormMatrices(*eigen_problem, eps, x, ctx));
800 :
801 688 : if (pc != jac)
802 : {
803 688 : LibmeshPetscCallQ(MatAssemblyBegin(jac, MAT_FINAL_ASSEMBLY));
804 688 : LibmeshPetscCallQ(MatAssemblyEnd(jac, MAT_FINAL_ASSEMBLY));
805 : }
806 688 : PetscFunctionReturn(PETSC_SUCCESS);
807 : }
808 :
809 3940 : LibmeshPetscCallQ(PetscObjectTypeCompare((PetscObject)pc, MATSHELL, &pisshell));
810 3940 : if ((jisshell || jismffd) && pisshell)
811 : {
812 : // Just assemble matrices and return
813 290 : LibmeshPetscCallQ(MatAssemblyBegin(jac, MAT_FINAL_ASSEMBLY));
814 290 : LibmeshPetscCallQ(MatAssemblyBegin(pc, MAT_FINAL_ASSEMBLY));
815 290 : LibmeshPetscCallQ(MatAssemblyEnd(jac, MAT_FINAL_ASSEMBLY));
816 290 : LibmeshPetscCallQ(MatAssemblyEnd(pc, MAT_FINAL_ASSEMBLY));
817 :
818 290 : PetscFunctionReturn(PETSC_SUCCESS);
819 : }
820 :
821 : // Jacobian and precond matrix are the same
822 3650 : if (jac == pc)
823 : {
824 240 : if (!pisshell)
825 240 : moosePetscSNESFormMatrixTag(
826 : snes, x, pc, sys.get_matrix_A(), ctx, eigen_nl.precondMatrixTag());
827 :
828 240 : PetscFunctionReturn(PETSC_SUCCESS);
829 : }
830 : else
831 : {
832 3410 : if (!jisshell && !jismffd && !pisshell) // We need to form both Jacobian and precond matrix
833 : {
834 0 : std::vector<Mat> mats = {jac, pc};
835 0 : std::vector<SparseMatrix<Number> *> libmesh_mats = {&sys.get_matrix_A(),
836 0 : &sys.get_precond_matrix()};
837 0 : std::set<TagID> tags = {eigen_nl.nonEigenMatrixTag(), eigen_nl.precondMatrixTag()};
838 0 : moosePetscSNESFormMatricesTags(snes, x, mats, libmesh_mats, ctx, tags);
839 0 : PetscFunctionReturn(PETSC_SUCCESS);
840 0 : }
841 3410 : if (!pisshell) // We need to form only precond matrix
842 : {
843 3410 : moosePetscSNESFormMatrixTag(
844 : snes, x, pc, sys.get_precond_matrix(), ctx, eigen_nl.precondMatrixTag());
845 3410 : LibmeshPetscCallQ(MatAssemblyBegin(jac, MAT_FINAL_ASSEMBLY));
846 3410 : LibmeshPetscCallQ(MatAssemblyEnd(jac, MAT_FINAL_ASSEMBLY));
847 3410 : PetscFunctionReturn(PETSC_SUCCESS);
848 : }
849 0 : if (!jisshell && !jismffd) // We need to form only Jacobian matrix
850 : {
851 0 : moosePetscSNESFormMatrixTag(
852 : snes, x, jac, sys.get_matrix_A(), ctx, eigen_nl.nonEigenMatrixTag());
853 0 : LibmeshPetscCallQ(MatAssemblyBegin(pc, MAT_FINAL_ASSEMBLY));
854 0 : LibmeshPetscCallQ(MatAssemblyEnd(pc, MAT_FINAL_ASSEMBLY));
855 0 : PetscFunctionReturn(PETSC_SUCCESS);
856 : }
857 : }
858 0 : PetscFunctionReturn(PETSC_SUCCESS);
859 : }
860 :
861 : PetscErrorCode
862 0 : mooseSlepcEigenFormJacobianB(SNES snes, Vec x, Mat jac, Mat pc, void * ctx)
863 : {
864 : PetscBool jshell, pshell;
865 : PetscBool jismffd;
866 :
867 : PetscFunctionBegin;
868 :
869 0 : EigenProblem * eigen_problem = static_cast<EigenProblem *>(ctx);
870 0 : NonlinearEigenSystem & eigen_nl = eigen_problem->getCurrentNonlinearEigenSystem();
871 0 : auto & sys = eigen_nl.sys();
872 :
873 : // If both jacobian and preconditioning are shell matrices,
874 : // and then assemble them and return
875 0 : LibmeshPetscCallQ(PetscObjectTypeCompare((PetscObject)jac, MATSHELL, &jshell));
876 0 : LibmeshPetscCallQ(PetscObjectTypeCompare((PetscObject)jac, MATMFFD, &jismffd));
877 0 : LibmeshPetscCallQ(PetscObjectTypeCompare((PetscObject)pc, MATSHELL, &pshell));
878 0 : if ((jshell || jismffd) && pshell)
879 : {
880 : // Just assemble matrices and return
881 0 : LibmeshPetscCallQ(MatAssemblyBegin(jac, MAT_FINAL_ASSEMBLY));
882 0 : LibmeshPetscCallQ(MatAssemblyBegin(pc, MAT_FINAL_ASSEMBLY));
883 0 : LibmeshPetscCallQ(MatAssemblyEnd(jac, MAT_FINAL_ASSEMBLY));
884 0 : LibmeshPetscCallQ(MatAssemblyEnd(pc, MAT_FINAL_ASSEMBLY));
885 :
886 0 : PetscFunctionReturn(PETSC_SUCCESS);
887 : }
888 :
889 0 : if (jac != pc && (!jshell && !jshell))
890 0 : SETERRQ(PetscObjectComm((PetscObject)snes),
891 : PETSC_ERR_ARG_INCOMP,
892 : "Jacobian and precond matrices should be the same for eigen kernels \n");
893 :
894 0 : moosePetscSNESFormMatrixTag(snes, x, pc, sys.get_matrix_B(), ctx, eigen_nl.eigenMatrixTag());
895 :
896 0 : if (eigen_problem->negativeSignEigenKernel())
897 : {
898 0 : LibmeshPetscCallQ(MatScale(pc, -1.));
899 : }
900 :
901 0 : PetscFunctionReturn(PETSC_SUCCESS);
902 : }
903 :
904 : void
905 17694 : moosePetscSNESFormFunction(SNES /*snes*/, Vec x, Vec r, void * ctx, TagID tag)
906 : {
907 17694 : EigenProblem * eigen_problem = static_cast<EigenProblem *>(ctx);
908 17694 : evaluateResidual(*eigen_problem, x, r, tag);
909 17694 : }
910 :
911 : PetscErrorCode
912 17210 : mooseSlepcEigenFormFunctionA(SNES snes, Vec x, Vec r, void * ctx)
913 : {
914 : PetscFunctionBegin;
915 :
916 17210 : EigenProblem * eigen_problem = static_cast<EigenProblem *>(ctx);
917 17210 : NonlinearEigenSystem & eigen_nl = eigen_problem->getCurrentNonlinearEigenSystem();
918 :
919 21140 : if (eigen_problem->solverParams(eigen_nl.number())._eigen_matrix_vector_mult &&
920 3930 : (eigen_problem->onLinearSolver() || eigen_problem->constantMatrices()))
921 : {
922 2810 : EPS eps = eigen_nl.getEPS();
923 : Mat A;
924 : ST st;
925 :
926 2810 : LibmeshPetscCallQ(mooseEPSFormMatrices(*eigen_problem, eps, x, ctx));
927 :
928 : // Rest ST state so that we can restrieve matrices
929 2810 : LibmeshPetscCallQ(EPSGetST(eps, &st));
930 2810 : LibmeshPetscCallQ(STResetMatrixState(st));
931 2810 : LibmeshPetscCallQ(EPSGetOperators(eps, &A, NULL));
932 :
933 2810 : LibmeshPetscCallQ(MatMult(A, x, r));
934 :
935 2810 : PetscFunctionReturn(PETSC_SUCCESS);
936 : }
937 :
938 14400 : moosePetscSNESFormFunction(snes, x, r, ctx, eigen_nl.nonEigenVectorTag());
939 :
940 14400 : PetscFunctionReturn(PETSC_SUCCESS);
941 : }
942 :
943 : PetscErrorCode
944 3558 : mooseSlepcEigenFormFunctionB(SNES snes, Vec x, Vec r, void * ctx)
945 : {
946 : PetscFunctionBegin;
947 :
948 3558 : EigenProblem * eigen_problem = static_cast<EigenProblem *>(ctx);
949 3558 : NonlinearEigenSystem & eigen_nl = eigen_problem->getCurrentNonlinearEigenSystem();
950 :
951 4662 : if (eigen_problem->solverParams(eigen_nl.number())._eigen_matrix_vector_mult &&
952 1104 : (eigen_problem->onLinearSolver() || eigen_problem->constantMatrices()))
953 : {
954 264 : EPS eps = eigen_nl.getEPS();
955 : ST st;
956 : Mat B;
957 :
958 264 : LibmeshPetscCallQ(mooseEPSFormMatrices(*eigen_problem, eps, x, ctx));
959 :
960 : // Rest ST state so that we can restrieve matrices
961 264 : LibmeshPetscCallQ(EPSGetST(eps, &st));
962 264 : LibmeshPetscCallQ(STResetMatrixState(st));
963 264 : LibmeshPetscCallQ(EPSGetOperators(eps, NULL, &B));
964 :
965 264 : LibmeshPetscCallQ(MatMult(B, x, r));
966 :
967 264 : if (eigen_problem->bxNormProvided())
968 : {
969 : // User has provided a postprocessor. We need it updated
970 0 : updateCurrentLocalSolution(eigen_nl.sys(), x);
971 0 : eigen_problem->execute(EXEC_LINEAR);
972 : }
973 : }
974 : else
975 3294 : moosePetscSNESFormFunction(snes, x, r, ctx, eigen_nl.eigenVectorTag());
976 :
977 3558 : if (eigen_problem->negativeSignEigenKernel())
978 : {
979 3558 : LibmeshPetscCallQ(VecScale(r, -1.));
980 : }
981 :
982 3558 : PetscFunctionReturn(PETSC_SUCCESS);
983 : }
984 :
985 : PetscErrorCode
986 15648 : mooseSlepcEigenFormFunctionAB(SNES /*snes*/, Vec x, Vec Ax, Vec Bx, void * ctx)
987 : {
988 : PetscFunctionBegin;
989 :
990 15648 : EigenProblem * eigen_problem = static_cast<EigenProblem *>(ctx);
991 15648 : NonlinearEigenSystem & eigen_nl = eigen_problem->getCurrentNonlinearEigenSystem();
992 15648 : auto & sys = eigen_nl.sys();
993 15648 : auto & dof_map = sys.get_dof_map();
994 :
995 18650 : if (eigen_problem->solverParams(eigen_nl.number())._eigen_matrix_vector_mult &&
996 3002 : (eigen_problem->onLinearSolver() || eigen_problem->constantMatrices()))
997 : {
998 2398 : EPS eps = eigen_nl.getEPS();
999 : ST st;
1000 : Mat A, B;
1001 :
1002 2398 : LibmeshPetscCallQ(mooseEPSFormMatrices(*eigen_problem, eps, x, ctx));
1003 :
1004 : // Rest ST state so that we can restrieve matrices
1005 2398 : LibmeshPetscCallQ(EPSGetST(eps, &st));
1006 2398 : LibmeshPetscCallQ(STResetMatrixState(st));
1007 :
1008 2398 : LibmeshPetscCallQ(EPSGetOperators(eps, &A, &B));
1009 :
1010 2398 : LibmeshPetscCallQ(MatMult(A, x, Ax));
1011 2398 : LibmeshPetscCallQ(MatMult(B, x, Bx));
1012 :
1013 2398 : if (eigen_problem->negativeSignEigenKernel())
1014 2398 : LibmeshPetscCallQ(VecScale(Bx, -1.));
1015 :
1016 2398 : if (eigen_problem->bxNormProvided())
1017 : {
1018 : // User has provided a postprocessor. We need it updated
1019 264 : updateCurrentLocalSolution(sys, x);
1020 264 : eigen_problem->execute(EXEC_LINEAR);
1021 : }
1022 :
1023 2398 : PetscFunctionReturn(PETSC_SUCCESS);
1024 : }
1025 :
1026 13250 : updateCurrentLocalSolution(sys, x);
1027 13250 : auto AX = createWrappedResidual(sys, Ax);
1028 13250 : auto BX = createWrappedResidual(sys, Bx);
1029 :
1030 26500 : eigen_problem->computeResidualAB(*sys.current_local_solution.get(),
1031 13250 : *AX,
1032 13250 : *BX,
1033 : eigen_nl.nonEigenVectorTag(),
1034 : eigen_nl.eigenVectorTag());
1035 :
1036 13250 : AX->close();
1037 13250 : BX->close();
1038 :
1039 13250 : if (dof_map.n_constrained_dofs())
1040 : {
1041 420 : PetscVector<Number> sub_Ax(Ax, sys.comm());
1042 420 : sys.copy_super_to_sub(*AX, sub_Ax);
1043 420 : PetscVector<Number> sub_Bx(Bx, sys.comm());
1044 420 : sys.copy_super_to_sub(*BX, sub_Bx);
1045 420 : }
1046 :
1047 13250 : if (eigen_problem->negativeSignEigenKernel())
1048 13250 : LibmeshPetscCallQ(VecScale(Bx, -1.));
1049 :
1050 13250 : PetscFunctionReturn(PETSC_SUCCESS);
1051 13250 : }
1052 :
1053 : PetscErrorCode
1054 2342 : mooseSlepcEigenFormNorm(SNES /*snes*/, Vec /*Bx*/, PetscReal * norm, void * ctx)
1055 : {
1056 : PetscFunctionBegin;
1057 2342 : auto * const eigen_problem = static_cast<EigenProblem *>(ctx);
1058 2342 : *norm = eigen_problem->formNorm();
1059 2342 : PetscFunctionReturn(PETSC_SUCCESS);
1060 : }
1061 :
1062 : void
1063 2044 : 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 2044 : LibmeshPetscCallA(
1077 : eigen_problem.comm().get(),
1078 : PetscObjectComposeFunction((PetscObject)mat,
1079 : "formJacobian",
1080 : eigen ? Moose::SlepcSupport::mooseSlepcEigenFormJacobianB
1081 : : Moose::SlepcSupport::mooseSlepcEigenFormJacobianA));
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 2044 : LibmeshPetscCallA(
1087 : eigen_problem.comm().get(),
1088 : PetscObjectComposeFunction((PetscObject)mat,
1089 : "formFunction",
1090 : eigen ? Moose::SlepcSupport::mooseSlepcEigenFormFunctionB
1091 : : Moose::SlepcSupport::mooseSlepcEigenFormFunctionA));
1092 :
1093 : // It's also beneficial to be able to evaluate both A and B residuals at once
1094 2044 : LibmeshPetscCallA(eigen_problem.comm().get(),
1095 : PetscObjectComposeFunction((PetscObject)mat,
1096 : "formFunctionAB",
1097 : Moose::SlepcSupport::mooseSlepcEigenFormFunctionAB));
1098 :
1099 : // Users may choose to provide a custom measure of the norm of B (Bx for a linear system)
1100 2044 : if (eigen_problem.bxNormProvided())
1101 84 : LibmeshPetscCallA(eigen_problem.comm().get(),
1102 : PetscObjectComposeFunction((PetscObject)mat,
1103 : "formNorm",
1104 : Moose::SlepcSupport::mooseSlepcEigenFormNorm));
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 2044 : LibmeshPetscCallA(eigen_problem.comm().get(),
1110 : PetscContainerCreate(eigen_problem.comm().get(), &container));
1111 2044 : LibmeshPetscCallA(eigen_problem.comm().get(),
1112 : PetscContainerSetPointer(container, &eigen_problem));
1113 2044 : LibmeshPetscCallA(
1114 : eigen_problem.comm().get(),
1115 : PetscObjectCompose((PetscObject)mat, "formJacobianCtx", (PetscObject)container));
1116 2044 : LibmeshPetscCallA(
1117 : eigen_problem.comm().get(),
1118 : PetscObjectCompose((PetscObject)mat, "formFunctionCtx", (PetscObject)container));
1119 2044 : if (eigen_problem.bxNormProvided())
1120 84 : LibmeshPetscCallA(eigen_problem.comm().get(),
1121 : PetscObjectCompose((PetscObject)mat, "formNormCtx", (PetscObject)container));
1122 :
1123 2044 : LibmeshPetscCallA(eigen_problem.comm().get(), PetscContainerDestroy(&container));
1124 2044 : }
1125 :
1126 : PetscErrorCode
1127 0 : mooseMatMult_Eigen(Mat mat, Vec x, Vec r)
1128 : {
1129 : PetscFunctionBegin;
1130 0 : void * ctx = nullptr;
1131 0 : LibmeshPetscCallQ(MatShellGetContext(mat, &ctx));
1132 :
1133 0 : if (!ctx)
1134 0 : mooseError("No context is set for shell matrix ");
1135 :
1136 0 : EigenProblem * eigen_problem = static_cast<EigenProblem *>(ctx);
1137 0 : NonlinearEigenSystem & eigen_nl = eigen_problem->getCurrentNonlinearEigenSystem();
1138 :
1139 0 : evaluateResidual(*eigen_problem, x, r, eigen_nl.eigenVectorTag());
1140 :
1141 0 : if (eigen_problem->negativeSignEigenKernel())
1142 0 : LibmeshPetscCallQ(VecScale(r, -1.));
1143 :
1144 0 : PetscFunctionReturn(PETSC_SUCCESS);
1145 : }
1146 :
1147 : PetscErrorCode
1148 0 : mooseMatMult_NonEigen(Mat mat, Vec x, Vec r)
1149 : {
1150 : PetscFunctionBegin;
1151 0 : void * ctx = nullptr;
1152 0 : LibmeshPetscCallQ(MatShellGetContext(mat, &ctx));
1153 :
1154 0 : if (!ctx)
1155 0 : mooseError("No context is set for shell matrix ");
1156 :
1157 0 : EigenProblem * eigen_problem = static_cast<EigenProblem *>(ctx);
1158 0 : NonlinearEigenSystem & eigen_nl = eigen_problem->getCurrentNonlinearEigenSystem();
1159 :
1160 0 : evaluateResidual(*eigen_problem, x, r, eigen_nl.nonEigenVectorTag());
1161 :
1162 0 : PetscFunctionReturn(PETSC_SUCCESS);
1163 : }
1164 :
1165 : void
1166 1168 : setOperationsForShellMat(EigenProblem & eigen_problem, Mat mat, bool eigen)
1167 : {
1168 1168 : LibmeshPetscCallA(eigen_problem.comm().get(), MatShellSetContext(mat, &eigen_problem));
1169 1168 : LibmeshPetscCallA(eigen_problem.comm().get(),
1170 : MatShellSetOperation(mat,
1171 : MATOP_MULT,
1172 : eigen ? (void (*)(void))mooseMatMult_Eigen
1173 : : (void (*)(void))mooseMatMult_NonEigen));
1174 1168 : }
1175 :
1176 : PETSC_EXTERN PetscErrorCode
1177 50 : registerPCToPETSc()
1178 : {
1179 : PetscFunctionBegin;
1180 :
1181 50 : LibmeshPetscCallQ(PCRegister("moosepc", PCCreate_MoosePC));
1182 :
1183 50 : PetscFunctionReturn(PETSC_SUCCESS);
1184 : }
1185 :
1186 : PETSC_EXTERN PetscErrorCode
1187 100 : PCCreate_MoosePC(PC pc)
1188 : {
1189 : PetscFunctionBegin;
1190 :
1191 100 : pc->ops->view = PCView_MoosePC;
1192 100 : pc->ops->destroy = PCDestroy_MoosePC;
1193 100 : pc->ops->setup = PCSetUp_MoosePC;
1194 100 : pc->ops->apply = PCApply_MoosePC;
1195 :
1196 100 : PetscFunctionReturn(PETSC_SUCCESS);
1197 : }
1198 :
1199 : PetscErrorCode
1200 100 : PCDestroy_MoosePC(PC /*pc*/)
1201 : {
1202 : PetscFunctionBegin;
1203 : /* We do not need to do anything right now, but later we may have some data we need to free here
1204 : */
1205 100 : PetscFunctionReturn(PETSC_SUCCESS);
1206 : }
1207 :
1208 : PetscErrorCode
1209 80 : PCView_MoosePC(PC /*pc*/, PetscViewer viewer)
1210 : {
1211 : PetscBool iascii;
1212 :
1213 : PetscFunctionBegin;
1214 80 : LibmeshPetscCallQ(PetscObjectTypeCompare((PetscObject)viewer, PETSCVIEWERASCII, &iascii));
1215 80 : if (iascii)
1216 80 : LibmeshPetscCallQ(PetscViewerASCIIPrintf(viewer, " %s\n", "moosepc"));
1217 :
1218 80 : PetscFunctionReturn(PETSC_SUCCESS);
1219 : }
1220 :
1221 : PetscErrorCode
1222 1230 : PCApply_MoosePC(PC pc, Vec x, Vec y)
1223 : {
1224 : void * ctx;
1225 : Mat Amat, Pmat;
1226 : PetscContainer container;
1227 :
1228 : PetscFunctionBegin;
1229 1230 : LibmeshPetscCallQ(PCGetOperators(pc, &Amat, &Pmat));
1230 1230 : LibmeshPetscCallQ(
1231 : PetscObjectQuery((PetscObject)Pmat, "formFunctionCtx", (PetscObject *)&container));
1232 1230 : if (container)
1233 1230 : LibmeshPetscCallQ(PetscContainerGetPointer(container, &ctx));
1234 : else
1235 0 : mooseError(" Can not find a context \n");
1236 :
1237 1230 : EigenProblem * eigen_problem = static_cast<EigenProblem *>(ctx);
1238 1230 : NonlinearEigenSystem & nl_eigen = eigen_problem->getCurrentNonlinearEigenSystem();
1239 1230 : auto preconditioner = nl_eigen.preconditioner();
1240 :
1241 1230 : if (!preconditioner)
1242 0 : mooseError("There is no moose preconditioner in nonlinear eigen system \n");
1243 :
1244 1230 : PetscVector<Number> x_vec(x, preconditioner->comm());
1245 1230 : PetscVector<Number> y_vec(y, preconditioner->comm());
1246 :
1247 1230 : preconditioner->apply(x_vec, y_vec);
1248 :
1249 1230 : PetscFunctionReturn(PETSC_SUCCESS);
1250 1230 : }
1251 :
1252 : PetscErrorCode
1253 360 : PCSetUp_MoosePC(PC pc)
1254 : {
1255 : void * ctx;
1256 : Mat Amat, Pmat;
1257 : PetscContainer container;
1258 :
1259 : PetscFunctionBegin;
1260 360 : LibmeshPetscCallQ(PCGetOperators(pc, &Amat, &Pmat));
1261 360 : LibmeshPetscCallQ(
1262 : PetscObjectQuery((PetscObject)Pmat, "formFunctionCtx", (PetscObject *)&container));
1263 360 : if (container)
1264 360 : LibmeshPetscCallQ(PetscContainerGetPointer(container, &ctx));
1265 : else
1266 0 : mooseError(" Can not find a context \n");
1267 :
1268 360 : EigenProblem * eigen_problem = static_cast<EigenProblem *>(ctx);
1269 360 : NonlinearEigenSystem & nl_eigen = eigen_problem->getCurrentNonlinearEigenSystem();
1270 360 : Preconditioner<Number> * preconditioner = nl_eigen.preconditioner();
1271 :
1272 360 : if (!preconditioner)
1273 0 : mooseError("There is no moose preconditioner in nonlinear eigen system \n");
1274 :
1275 360 : if (!preconditioner->initialized())
1276 50 : preconditioner->init();
1277 :
1278 360 : preconditioner->setup();
1279 :
1280 360 : PetscFunctionReturn(PETSC_SUCCESS);
1281 : }
1282 :
1283 : PetscErrorCode
1284 2988 : mooseSlepcStoppingTest(EPS eps,
1285 : PetscInt its,
1286 : PetscInt max_it,
1287 : PetscInt nconv,
1288 : PetscInt nev,
1289 : EPSConvergedReason * reason,
1290 : void * ctx)
1291 : {
1292 2988 : EigenProblem * eigen_problem = static_cast<EigenProblem *>(ctx);
1293 :
1294 : PetscFunctionBegin;
1295 2988 : 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 2988 : if (eigen_problem->doFreePowerIteration() && its == max_it && *reason <= 0)
1304 : {
1305 534 : *reason = EPS_CONVERGED_USER;
1306 534 : eps->nconv = 1;
1307 : }
1308 2988 : PetscFunctionReturn(PETSC_SUCCESS);
1309 : }
1310 :
1311 : PetscErrorCode
1312 8854 : mooseSlepcEPSGetSNES(EPS eps, SNES * snes)
1313 : {
1314 : PetscBool same, nonlinear;
1315 :
1316 : PetscFunctionBegin;
1317 8854 : LibmeshPetscCallQ(PetscObjectTypeCompare((PetscObject)eps, EPSPOWER, &same));
1318 :
1319 8854 : if (!same)
1320 0 : mooseError("It is not eps power, and there is no snes");
1321 :
1322 8854 : LibmeshPetscCallQ(EPSPowerGetNonlinear(eps, &nonlinear));
1323 :
1324 8854 : if (!nonlinear)
1325 0 : mooseError("It is not a nonlinear eigen solver");
1326 :
1327 8854 : LibmeshPetscCallQ(EPSPowerGetSNES(eps, snes));
1328 :
1329 8854 : PetscFunctionReturn(PETSC_SUCCESS);
1330 : }
1331 :
1332 : PetscErrorCode
1333 1246 : mooseSlepcEPSSNESSetUpOptionPrefix(EPS eps)
1334 : {
1335 : SNES snes;
1336 1246 : const char * prefix = nullptr;
1337 :
1338 : PetscFunctionBegin;
1339 1246 : LibmeshPetscCallQ(mooseSlepcEPSGetSNES(eps, &snes));
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 1246 : LibmeshPetscCallQ(PetscObjectGetOptionsPrefix((PetscObject)eps, &prefix));
1344 : // Set option prefix to SNES
1345 1246 : LibmeshPetscCallQ(SNESSetOptionsPrefix(snes, prefix));
1346 :
1347 1246 : PetscFunctionReturn(PETSC_SUCCESS);
1348 : }
1349 :
1350 : PetscErrorCode
1351 100 : mooseSlepcEPSSNESSetCustomizePC(EPS eps)
1352 : {
1353 : SNES snes;
1354 : KSP ksp;
1355 : PC pc;
1356 :
1357 : PetscFunctionBegin;
1358 : // Get SNES from EPS
1359 100 : LibmeshPetscCallQ(mooseSlepcEPSGetSNES(eps, &snes));
1360 : // Get KSP from SNES
1361 100 : LibmeshPetscCallQ(SNESGetKSP(snes, &ksp));
1362 : // Get PC from KSP
1363 100 : LibmeshPetscCallQ(KSPGetPC(ksp, &pc));
1364 : // Set PC type
1365 100 : LibmeshPetscCallQ(PCSetType(pc, "moosepc"));
1366 100 : PetscFunctionReturn(PETSC_SUCCESS);
1367 : }
1368 :
1369 : PetscErrorCode
1370 1246 : mooseSlepcEPSSNESKSPSetPCSide(FEProblemBase & problem, EPS eps)
1371 : {
1372 : SNES snes;
1373 : KSP ksp;
1374 :
1375 : PetscFunctionBegin;
1376 : // Get SNES from EPS
1377 1246 : LibmeshPetscCallQ(mooseSlepcEPSGetSNES(eps, &snes));
1378 : // Get KSP from SNES
1379 1246 : LibmeshPetscCallQ(SNESGetKSP(snes, &ksp));
1380 :
1381 1246 : Moose::PetscSupport::petscSetDefaultPCSide(problem, ksp);
1382 :
1383 1246 : Moose::PetscSupport::petscSetDefaultKSPNormType(problem, ksp);
1384 1246 : PetscFunctionReturn(PETSC_SUCCESS);
1385 : }
1386 :
1387 : PetscErrorCode
1388 5112 : mooseSlepcEPSMonitor(EPS eps,
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 :
1400 : PetscFunctionBegin;
1401 5112 : EigenProblem * eigen_problem = static_cast<EigenProblem *>(mctx);
1402 5112 : auto & console = eigen_problem->console();
1403 :
1404 5112 : auto inverse = eigen_problem->outputInverseEigenvalue();
1405 5112 : LibmeshPetscCallQ(EPSGetST(eps, &st));
1406 5112 : eigenr = eigr[0];
1407 5112 : eigeni = eigi[0];
1408 : // Make the eigenvalue consistent with shift type
1409 5112 : LibmeshPetscCallQ(STBackTransform(st, 1, &eigenr, &eigeni));
1410 :
1411 5112 : auto eigenvalue = inverse ? 1.0 / eigenr : eigenr;
1412 :
1413 : // The term "k-eigenvalue" is adopted from the neutronics community.
1414 5112 : console << " Iteration " << its << std::setprecision(10) << std::fixed
1415 5112 : << (inverse ? " k-eigenvalue = " : " eigenvalue = ") << eigenvalue << std::endl;
1416 :
1417 5112 : PetscFunctionReturn(PETSC_SUCCESS);
1418 : }
1419 :
1420 : } // namespace SlepcSupport
1421 : } // namespace moose
1422 :
1423 : #endif // LIBMESH_HAVE_SLEPC
|
