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