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 : #include "EigenProblem.h"
13 : #include "Assembly.h"
14 : #include "AuxiliarySystem.h"
15 : #include "DisplacedProblem.h"
16 : #include "NonlinearEigenSystem.h"
17 : #include "SlepcSupport.h"
18 : #include "RandomData.h"
19 : #include "OutputWarehouse.h"
20 : #include "Function.h"
21 : #include "MooseVariableScalar.h"
22 : #include "UserObject.h"
23 :
24 : // libMesh includes
25 : #include "libmesh/system.h"
26 : #include "libmesh/eigen_solver.h"
27 : #include "libmesh/enum_eigen_solver_type.h"
28 :
29 : // Needed for LIBMESH_CHECK_ERR
30 : using libMesh::PetscSolverException;
31 :
32 : registerMooseObject("MooseApp", EigenProblem);
33 :
34 : InputParameters
35 15389 : EigenProblem::validParams()
36 : {
37 15389 : InputParameters params = FEProblemBase::validParams();
38 15389 : params.addClassDescription("Problem object for solving an eigenvalue problem.");
39 46167 : params.addParam<bool>("negative_sign_eigen_kernel",
40 30778 : true,
41 : "Whether or not to use a negative sign for eigenvalue kernels. "
42 : "Using a negative sign makes eigenvalue kernels consistent with "
43 : "a nonlinear solver");
44 :
45 46167 : params.addParam<unsigned int>(
46 : "active_eigen_index",
47 30778 : 0,
48 : "Which eigenvector is used to compute residual and also associated to nonlinear variable");
49 15389 : params.addParam<PostprocessorName>("bx_norm", "A postprocessor describing the norm of Bx");
50 :
51 15389 : params.addParamNamesToGroup("negative_sign_eigen_kernel active_eigen_index bx_norm",
52 : "Eigenvalue solve");
53 :
54 15389 : return params;
55 0 : }
56 :
57 562 : EigenProblem::EigenProblem(const InputParameters & parameters)
58 : : FEProblemBase(parameters)
59 : #ifdef LIBMESH_HAVE_SLEPC
60 : ,
61 : // By default, we want to compute an eigenvalue only (smallest or largest)
62 562 : _n_eigen_pairs_required(1),
63 562 : _generalized_eigenvalue_problem(false),
64 562 : _negative_sign_eigen_kernel(getParam<bool>("negative_sign_eigen_kernel")),
65 562 : _active_eigen_index(getParam<unsigned int>("active_eigen_index")),
66 562 : _do_free_power_iteration(false),
67 562 : _output_inverse_eigenvalue(false),
68 562 : _on_linear_solver(false),
69 562 : _matrices_formed(false),
70 562 : _constant_matrices(false),
71 562 : _has_normalization(false),
72 562 : _normal_factor(1.0),
73 562 : _first_solve(declareRestartableData<bool>("first_solve", true)),
74 1124 : _bx_norm_name(isParamValid("bx_norm")
75 562 : ? std::make_optional(getParam<PostprocessorName>("bx_norm"))
76 1124 : : std::nullopt)
77 : #endif
78 : {
79 : #ifdef LIBMESH_HAVE_SLEPC
80 562 : if (_nl_sys_names.size() > 1)
81 0 : paramError("nl_sys_names",
82 : "eigen problems do not currently support multiple nonlinear eigen systems");
83 562 : if (_linear_sys_names.size())
84 0 : paramError("linear_sys_names", "EigenProblem only works with a single nonlinear eigen system");
85 :
86 1124 : for (const auto i : index_range(_nl_sys_names))
87 : {
88 562 : const auto & sys_name = _nl_sys_names[i];
89 562 : auto & nl = _nl[i];
90 562 : nl = std::make_shared<NonlinearEigenSystem>(*this, sys_name);
91 562 : _nl_eigen = std::dynamic_pointer_cast<NonlinearEigenSystem>(nl);
92 562 : _current_nl_sys = nl.get();
93 562 : _solver_systems[i] = std::dynamic_pointer_cast<SolverSystem>(nl);
94 562 : nl->system().prefer_hash_table_matrix_assembly(_use_hash_table_matrix_assembly);
95 : }
96 :
97 562 : _aux = std::make_shared<AuxiliarySystem>(*this, "aux0");
98 :
99 562 : newAssemblyArray(_solver_systems);
100 :
101 562 : FEProblemBase::initNullSpaceVectors(parameters, _nl);
102 :
103 562 : es().parameters.set<EigenProblem *>("_eigen_problem") = this;
104 : #else
105 : mooseError("Need to install SLEPc to solve eigenvalue problems, please reconfigure\n");
106 : #endif /* LIBMESH_HAVE_SLEPC */
107 :
108 : // SLEPc older than 3.13.0 can not take initial guess from moose
109 : #if PETSC_RELEASE_LESS_THAN(3, 13, 0)
110 : mooseDeprecated(
111 : "Please use SLEPc-3.13.0 or higher. Old versions of SLEPc likely produce bad convergence");
112 : #endif
113 : // Create extra vectors and matrices if any
114 562 : createTagVectors();
115 :
116 : // Create extra solution vectors if any
117 562 : createTagSolutions();
118 562 : }
119 :
120 : #ifdef LIBMESH_HAVE_SLEPC
121 : void
122 558 : EigenProblem::setEigenproblemType(Moose::EigenProblemType eigen_problem_type)
123 : {
124 558 : switch (eigen_problem_type)
125 : {
126 28 : case Moose::EPT_HERMITIAN:
127 28 : _nl_eigen->sys().set_eigenproblem_type(libMesh::HEP);
128 28 : _generalized_eigenvalue_problem = false;
129 28 : break;
130 :
131 12 : case Moose::EPT_NON_HERMITIAN:
132 12 : _nl_eigen->sys().set_eigenproblem_type(libMesh::NHEP);
133 12 : _generalized_eigenvalue_problem = false;
134 12 : break;
135 :
136 0 : case Moose::EPT_GEN_HERMITIAN:
137 0 : _nl_eigen->sys().set_eigenproblem_type(libMesh::GHEP);
138 0 : _generalized_eigenvalue_problem = true;
139 0 : break;
140 :
141 0 : case Moose::EPT_GEN_INDEFINITE:
142 0 : _nl_eigen->sys().set_eigenproblem_type(libMesh::GHIEP);
143 0 : _generalized_eigenvalue_problem = true;
144 0 : break;
145 :
146 518 : case Moose::EPT_GEN_NON_HERMITIAN:
147 518 : _nl_eigen->sys().set_eigenproblem_type(libMesh::GNHEP);
148 518 : _generalized_eigenvalue_problem = true;
149 518 : break;
150 :
151 0 : case Moose::EPT_POS_GEN_NON_HERMITIAN:
152 0 : mooseError("libMesh does not support EPT_POS_GEN_NON_HERMITIAN currently \n");
153 : break;
154 :
155 0 : case Moose::EPT_SLEPC_DEFAULT:
156 0 : _generalized_eigenvalue_problem = false;
157 0 : break;
158 :
159 0 : default:
160 0 : mooseError("Unknown eigen solver type \n");
161 : }
162 558 : }
163 :
164 : void
165 4796 : EigenProblem::execute(const ExecFlagType & exec_type)
166 : {
167 4796 : if (exec_type == EXEC_INITIAL && !_app.isRestarting())
168 : // we need to scale the solution properly and we can do this only all initial setup of
169 : // depending objects by the residual evaluations has been done to this point.
170 506 : preScaleEigenVector(std::pair<Real, Real>(_initial_eigenvalue, 0));
171 :
172 4796 : FEProblemBase::execute(exec_type);
173 4796 : }
174 :
175 : void
176 3481 : EigenProblem::computeJacobianTag(const NumericVector<Number> & soln,
177 : SparseMatrix<Number> & jacobian,
178 : TagID tag)
179 : {
180 3481 : TIME_SECTION("computeJacobianTag", 3);
181 :
182 : // Disassociate the default tags because we will associate vectors with only the
183 : // specific system tags that we need for this instance
184 3481 : _nl_eigen->disassociateDefaultMatrixTags();
185 :
186 : // Clear FE tags and first add the specific tag associated with the Jacobian
187 3481 : _fe_matrix_tags.clear();
188 3481 : _fe_matrix_tags.insert(tag);
189 :
190 : // Add any other user-added matrix tags if they have associated matrices
191 3481 : const auto & matrix_tags = getMatrixTags();
192 21328 : for (const auto & matrix_tag : matrix_tags)
193 17847 : if (_nl_eigen->hasMatrix(matrix_tag.second))
194 442 : _fe_matrix_tags.insert(matrix_tag.second);
195 :
196 3481 : _nl_eigen->setSolution(soln);
197 :
198 3481 : _nl_eigen->associateMatrixToTag(jacobian, tag);
199 :
200 3481 : setCurrentNonlinearSystem(_nl_eigen->number());
201 3481 : computeJacobianTags(_fe_matrix_tags);
202 :
203 3481 : _nl_eigen->disassociateMatrixFromTag(jacobian, tag);
204 3481 : }
205 :
206 : void
207 350 : EigenProblem::computeMatricesTags(const NumericVector<Number> & soln,
208 : const std::vector<SparseMatrix<Number> *> & jacobians,
209 : const std::set<TagID> & tags)
210 : {
211 350 : TIME_SECTION("computeMatricesTags", 3);
212 :
213 350 : if (jacobians.size() != tags.size())
214 0 : mooseError("The number of matrices ",
215 0 : jacobians.size(),
216 : " does not equal the number of tags ",
217 0 : tags.size());
218 :
219 : // Disassociate the default tags because we will associate vectors with only the
220 : // specific system tags that we need for this instance
221 350 : _nl_eigen->disassociateDefaultMatrixTags();
222 :
223 350 : _fe_matrix_tags.clear();
224 :
225 350 : _nl_eigen->setSolution(soln);
226 :
227 350 : unsigned int i = 0;
228 1050 : for (auto tag : tags)
229 700 : _nl_eigen->associateMatrixToTag(*(jacobians[i++]), tag);
230 :
231 350 : setCurrentNonlinearSystem(_nl_eigen->number());
232 350 : computeJacobianTags(tags);
233 :
234 350 : i = 0;
235 1050 : for (auto tag : tags)
236 700 : _nl_eigen->disassociateMatrixFromTag(*(jacobians[i++]), tag);
237 350 : }
238 :
239 : void
240 360 : EigenProblem::computeJacobianBlocks(std::vector<JacobianBlock *> & blocks,
241 : const unsigned int nl_sys_num)
242 : {
243 360 : TIME_SECTION("computeJacobianBlocks", 3);
244 360 : setCurrentNonlinearSystem(nl_sys_num);
245 :
246 360 : if (_displaced_problem)
247 0 : computeSystems(EXEC_PRE_DISPLACE);
248 :
249 360 : computeSystems(EXEC_NONLINEAR);
250 :
251 360 : _currently_computing_jacobian = true;
252 :
253 360 : _current_nl_sys->computeJacobianBlocks(blocks, {_nl_eigen->precondMatrixTag()});
254 :
255 360 : _currently_computing_jacobian = false;
256 360 : }
257 :
258 : void
259 22 : EigenProblem::computeJacobianAB(const NumericVector<Number> & soln,
260 : SparseMatrix<Number> & jacobianA,
261 : SparseMatrix<Number> & jacobianB,
262 : TagID tagA,
263 : TagID tagB)
264 : {
265 22 : TIME_SECTION("computeJacobianAB", 3);
266 :
267 : // Disassociate the default tags because we will associate vectors with only the
268 : // specific system tags that we need for this instance
269 22 : _nl_eigen->disassociateDefaultMatrixTags();
270 :
271 : // Clear FE tags and first add the specific tags associated with the Jacobian
272 22 : _fe_matrix_tags.clear();
273 22 : _fe_matrix_tags.insert(tagA);
274 22 : _fe_matrix_tags.insert(tagB);
275 :
276 : // Add any other user-added matrix tags if they have associated matrices
277 22 : const auto & matrix_tags = getMatrixTags();
278 132 : for (const auto & matrix_tag : matrix_tags)
279 110 : if (_nl_eigen->hasMatrix(matrix_tag.second))
280 0 : _fe_matrix_tags.insert(matrix_tag.second);
281 :
282 22 : _nl_eigen->setSolution(soln);
283 :
284 22 : _nl_eigen->associateMatrixToTag(jacobianA, tagA);
285 22 : _nl_eigen->associateMatrixToTag(jacobianB, tagB);
286 :
287 22 : setCurrentNonlinearSystem(_nl_eigen->number());
288 22 : computeJacobianTags(_fe_matrix_tags);
289 :
290 22 : _nl_eigen->disassociateMatrixFromTag(jacobianA, tagA);
291 22 : _nl_eigen->disassociateMatrixFromTag(jacobianB, tagB);
292 22 : }
293 :
294 : void
295 17487 : EigenProblem::computeResidualTag(const NumericVector<Number> & soln,
296 : NumericVector<Number> & residual,
297 : TagID tag)
298 : {
299 17487 : TIME_SECTION("computeResidualTag", 3);
300 :
301 : // Disassociate the default tags because we will associate vectors with only the
302 : // specific system tags that we need for this instance
303 17487 : _nl_eigen->disassociateDefaultVectorTags();
304 :
305 : // add the specific tag associated with the residual
306 : mooseAssert(_fe_vector_tags.empty(), "This should be empty indicating a clean starting state");
307 17487 : _fe_vector_tags.insert(tag);
308 :
309 : // Add any other user-added vector residual tags if they have associated vectors
310 17487 : const auto & residual_vector_tags = getVectorTags(Moose::VECTOR_TAG_RESIDUAL);
311 87765 : for (const auto & vector_tag : residual_vector_tags)
312 70278 : if (_nl_eigen->hasVector(vector_tag._id))
313 330 : _fe_vector_tags.insert(vector_tag._id);
314 :
315 17487 : _nl_eigen->associateVectorToTag(residual, tag);
316 :
317 17487 : _nl_eigen->setSolution(soln);
318 :
319 17487 : setCurrentNonlinearSystem(_nl_eigen->number());
320 17487 : computeResidualTags(_fe_vector_tags);
321 17487 : _fe_vector_tags.clear();
322 :
323 17487 : _nl_eigen->disassociateVectorFromTag(residual, tag);
324 17487 : }
325 :
326 : void
327 12889 : EigenProblem::computeResidualAB(const NumericVector<Number> & soln,
328 : NumericVector<Number> & residualA,
329 : NumericVector<Number> & residualB,
330 : TagID tagA,
331 : TagID tagB)
332 : {
333 12889 : TIME_SECTION("computeResidualAB", 3);
334 :
335 : // Disassociate the default tags because we will associate vectors with only the
336 : // specific system tags that we need for this instance
337 12889 : _nl_eigen->disassociateDefaultVectorTags();
338 :
339 : // add the specific tags associated with the residual
340 : mooseAssert(_fe_vector_tags.empty(), "This should be empty indicating a clean starting state");
341 12889 : _fe_vector_tags.insert(tagA);
342 12889 : _fe_vector_tags.insert(tagB);
343 :
344 : // Add any other user-added vector residual tags if they have associated vectors
345 12889 : const auto & residual_vector_tags = getVectorTags(Moose::VECTOR_TAG_RESIDUAL);
346 64821 : for (const auto & vector_tag : residual_vector_tags)
347 51932 : if (_nl_eigen->hasVector(vector_tag._id))
348 376 : _fe_vector_tags.insert(vector_tag._id);
349 :
350 12889 : _nl_eigen->associateVectorToTag(residualA, tagA);
351 12889 : _nl_eigen->associateVectorToTag(residualB, tagB);
352 :
353 12889 : _nl_eigen->setSolution(soln);
354 :
355 12889 : computeResidualTags(_fe_vector_tags);
356 12889 : _fe_vector_tags.clear();
357 :
358 12889 : _nl_eigen->disassociateVectorFromTag(residualA, tagA);
359 12889 : _nl_eigen->disassociateVectorFromTag(residualB, tagB);
360 12889 : }
361 :
362 : Real
363 132 : EigenProblem::computeResidualL2Norm()
364 : {
365 528 : computeResidualAB(*_nl_eigen->currentSolution(),
366 132 : _nl_eigen->residualVectorAX(),
367 132 : _nl_eigen->residualVectorBX(),
368 132 : _nl_eigen->nonEigenVectorTag(),
369 132 : _nl_eigen->eigenVectorTag());
370 :
371 132 : Real eigenvalue = 1.0;
372 :
373 132 : if (_active_eigen_index < _nl_eigen->getNumConvergedEigenvalues())
374 110 : eigenvalue = _nl_eigen->getConvergedEigenvalue(_active_eigen_index).first;
375 :
376 : // Scale BX with eigenvalue
377 132 : _nl_eigen->residualVectorBX() *= eigenvalue;
378 :
379 : // Compute entire residual
380 132 : if (_negative_sign_eigen_kernel)
381 132 : _nl_eigen->residualVectorAX() += _nl_eigen->residualVectorBX();
382 : else
383 0 : _nl_eigen->residualVectorAX() -= _nl_eigen->residualVectorBX();
384 :
385 132 : return _nl_eigen->residualVectorAX().l2_norm();
386 : }
387 :
388 : void
389 1048 : EigenProblem::adjustEigenVector(const Real value, bool scaling)
390 : {
391 1048 : std::vector<VariableName> var_names = getVariableNames();
392 3109 : for (auto & vn : var_names)
393 : {
394 2061 : MooseVariableBase * var = nullptr;
395 2061 : if (hasScalarVariable(vn))
396 46 : var = &getScalarVariable(0, vn);
397 : else
398 2015 : var = &getVariable(0, vn);
399 : // Do operations for only eigen variable
400 2061 : if (var->eigen())
401 2674 : for (unsigned int vc = 0; vc < var->count(); ++vc)
402 : {
403 1446 : std::set<dof_id_type> var_indices;
404 1446 : _nl_eigen->system().local_dof_indices(var->number() + vc, var_indices);
405 92336 : for (const auto & dof : var_indices)
406 90890 : _nl_eigen->solution().set(dof, scaling ? (_nl_eigen->solution()(dof) * value) : value);
407 1446 : }
408 : }
409 :
410 1048 : _nl_eigen->solution().close();
411 1048 : _nl_eigen->update();
412 1048 : }
413 :
414 : void
415 516 : EigenProblem::scaleEigenvector(const Real scaling_factor)
416 : {
417 516 : adjustEigenVector(scaling_factor, true);
418 516 : }
419 :
420 : void
421 532 : EigenProblem::initEigenvector(const Real initial_value)
422 : {
423 : // Yaqi's note: the following code will set a flat solution for lagrange and
424 : // constant monomial variables. For the first or higher order elemental variables,
425 : // the solution is not flat. Fortunately, the initial guess does not affect
426 : // the final solution as long as it is not perpendicular to the true solution.
427 : // We, in general, do not need to worry about that.
428 :
429 532 : adjustEigenVector(initial_value, false);
430 532 : }
431 :
432 : void
433 704 : EigenProblem::preScaleEigenVector(const std::pair<Real, Real> & eig)
434 : {
435 : // pre-scale the solution to make sure ||Bx||_2 is equal to inverse of eigenvalue
436 1408 : computeResidualTag(
437 2112 : *_nl_eigen->currentSolution(), _nl_eigen->residualVectorBX(), _nl_eigen->eigenVectorTag());
438 :
439 : // Eigenvalue magnitude
440 704 : Real v = std::sqrt(eig.first * eig.first + eig.second * eig.second);
441 : // Scaling factor
442 704 : Real factor = 1 / v / (bxNormProvided() ? formNorm() : _nl_eigen->residualVectorBX().l2_norm());
443 : // Scale eigenvector
444 704 : if (!MooseUtils::absoluteFuzzyEqual(factor, 1))
445 506 : scaleEigenvector(factor);
446 704 : }
447 :
448 : void
449 704 : EigenProblem::postScaleEigenVector()
450 : {
451 704 : if (_has_normalization)
452 : {
453 : Real v;
454 10 : if (_normal_factor == std::numeric_limits<Real>::max())
455 : {
456 0 : if (_active_eigen_index >= _nl_eigen->getNumConvergedEigenvalues())
457 0 : mooseError("Number of converged eigenvalues ",
458 0 : _nl_eigen->getNumConvergedEigenvalues(),
459 : " but you required eigenvalue ",
460 0 : _active_eigen_index);
461 :
462 : // when normal factor is not provided, we use the inverse of the norm of
463 : // the active eigenvalue for normalization
464 0 : auto eig = _nl_eigen->getAllConvergedEigenvalues()[_active_eigen_index];
465 0 : v = 1 / std::sqrt(eig.first * eig.first + eig.second * eig.second);
466 : }
467 : else
468 10 : v = _normal_factor;
469 :
470 10 : Real c = getPostprocessorValueByName(_normalization);
471 :
472 : // We scale SLEPc eigen vector here, so we need to scale it back for optimal
473 : // convergence if we call EPS solver again
474 : mooseAssert(v != 0., "normal factor can not be zero");
475 :
476 10 : unsigned int itr = 0;
477 :
478 20 : while (!MooseUtils::relativeFuzzyEqual(v, c))
479 : {
480 : // If postprocessor is not defined on eigen variables, scaling might not work
481 10 : if (itr > 10)
482 0 : mooseError("Can not scale eigenvector to the required factor ",
483 : v,
484 : " please check if postprocessor is defined on only eigen variables");
485 :
486 : mooseAssert(c != 0., "postprocessor value used for scaling can not be zero");
487 :
488 10 : scaleEigenvector(v / c);
489 :
490 : // update all aux variables and user objects on linear
491 10 : execute(EXEC_LINEAR);
492 :
493 10 : c = getPostprocessorValueByName(_normalization);
494 :
495 10 : itr++;
496 : }
497 : }
498 704 : }
499 :
500 : void
501 554 : EigenProblem::checkProblemIntegrity()
502 : {
503 554 : FEProblemBase::checkProblemIntegrity();
504 554 : _nl_eigen->checkIntegrity();
505 546 : if (_bx_norm_name)
506 : {
507 36 : if (!isNonlinearEigenvalueSolver(0))
508 0 : paramWarning("bx_norm", "This parameter is only used for nonlinear solve types");
509 36 : else if (auto & pp = getUserObjectBase(_bx_norm_name.value());
510 36 : !pp.getExecuteOnEnum().contains(EXEC_LINEAR))
511 4 : pp.paramError("execute_on",
512 : "If providing the Bx norm, this postprocessor must execute on linear e.g. "
513 : "during residual evaluations");
514 : }
515 542 : }
516 :
517 : void
518 529 : EigenProblem::doFreeNonlinearPowerIterations(unsigned int free_power_iterations)
519 : {
520 : mooseAssert(_current_nl_sys, "This needs to be non-null");
521 :
522 529 : doFreePowerIteration(true);
523 : // Set free power iterations
524 529 : Moose::SlepcSupport::setFreeNonlinearPowerIterations(free_power_iterations);
525 :
526 : // Call solver
527 529 : _current_nl_sys->solve();
528 529 : _current_nl_sys->update();
529 :
530 : // Clear free power iterations
531 529 : auto executioner = getMooseApp().getExecutioner();
532 529 : if (executioner)
533 529 : Moose::SlepcSupport::clearFreeNonlinearPowerIterations(executioner->parameters());
534 : else
535 0 : mooseError("There is no executioner for this moose app");
536 :
537 529 : doFreePowerIteration(false);
538 529 : }
539 :
540 : void
541 714 : EigenProblem::solve(const unsigned int nl_sys_num)
542 : {
543 : #if !PETSC_RELEASE_LESS_THAN(3, 12, 0)
544 : // Master has the default database
545 714 : if (!_app.isUltimateMaster())
546 132 : LibmeshPetscCall(PetscOptionsPush(_petsc_option_data_base));
547 : #endif
548 :
549 714 : setCurrentNonlinearSystem(nl_sys_num);
550 :
551 714 : if (_solve)
552 : {
553 704 : TIME_SECTION("solve", 1);
554 :
555 : // Set necessary slepc callbacks
556 : // We delay this call as much as possible because libmesh
557 : // could rebuild matrices due to mesh changes or something else.
558 704 : _nl_eigen->attachSLEPcCallbacks();
559 :
560 : // If there is an eigenvalue, we scale 1/|Bx| to eigenvalue
561 704 : if (_active_eigen_index < _nl_eigen->getNumConvergedEigenvalues())
562 : {
563 198 : std::pair<Real, Real> eig = _nl_eigen->getConvergedEigenvalue(_active_eigen_index);
564 198 : preScaleEigenVector(eig);
565 : }
566 :
567 1353 : if (isNonlinearEigenvalueSolver(nl_sys_num) &&
568 649 : solverParams(nl_sys_num)._eigen_solve_type != Moose::EST_NONLINEAR_POWER)
569 : {
570 : // Let do an initial solve if a nonlinear eigen solver but not power is used.
571 : // The initial solver is a Inverse Power, and it is used to compute a good initial
572 : // guess for Newton
573 638 : if (solverParams(nl_sys_num)._free_power_iterations && _first_solve)
574 : {
575 509 : _console << std::endl << " -------------------------------" << std::endl;
576 509 : _console << " Free power iteration starts ..." << std::endl;
577 509 : _console << " -------------------------------" << std::endl << std::endl;
578 509 : doFreeNonlinearPowerIterations(solverParams(nl_sys_num)._free_power_iterations);
579 509 : _first_solve = false;
580 : }
581 :
582 : // Let us do extra power iterations here if necessary
583 638 : if (solverParams(nl_sys_num)._extra_power_iterations)
584 : {
585 20 : _console << std::endl << " --------------------------------------" << std::endl;
586 20 : _console << " Extra Free power iteration starts ..." << std::endl;
587 20 : _console << " --------------------------------------" << std::endl << std::endl;
588 20 : doFreeNonlinearPowerIterations(solverParams(nl_sys_num)._extra_power_iterations);
589 : }
590 : }
591 :
592 : // We print this for only nonlinear solver
593 704 : if (isNonlinearEigenvalueSolver(nl_sys_num))
594 : {
595 649 : _console << std::endl << " -------------------------------------" << std::endl;
596 :
597 649 : if (solverParams(nl_sys_num)._eigen_solve_type != Moose::EST_NONLINEAR_POWER)
598 638 : _console << " Nonlinear Newton iteration starts ..." << std::endl;
599 : else
600 11 : _console << " Nonlinear power iteration starts ..." << std::endl;
601 :
602 649 : _console << " -------------------------------------" << std::endl << std::endl;
603 : }
604 :
605 704 : _current_nl_sys->solve();
606 704 : _current_nl_sys->update();
607 :
608 : // with PJFNKMO solve type, we need to evaluate the linear objects to bring them up-to-date
609 704 : if (solverParams(nl_sys_num)._eigen_solve_type == Moose::EST_PJFNKMO)
610 86 : execute(EXEC_LINEAR);
611 :
612 : // Scale eigen vector if users ask
613 704 : postScaleEigenVector();
614 704 : }
615 :
616 : #if !PETSC_RELEASE_LESS_THAN(3, 12, 0)
617 714 : if (!_app.isUltimateMaster())
618 132 : LibmeshPetscCall(PetscOptionsPop());
619 : #endif
620 :
621 : // sync solutions in displaced problem
622 714 : if (_displaced_problem)
623 33 : _displaced_problem->syncSolutions();
624 :
625 : // Reset the matrix flag, so that we reform matrix in next picard iteration
626 714 : _matrices_formed = false;
627 714 : }
628 :
629 : void
630 10 : EigenProblem::setNormalization(const PostprocessorName & pp, const Real value)
631 : {
632 10 : _has_normalization = true;
633 10 : _normalization = pp;
634 10 : _normal_factor = value;
635 10 : }
636 :
637 : void
638 554 : EigenProblem::init()
639 : {
640 : #if PETSC_RELEASE_LESS_THAN(3, 13, 0)
641 : // Prior to Slepc 3.13 we did not have a nonlinear eigenvalue solver so we must always assemble
642 : // before the solve
643 : _nl_eigen->sys().attach_assemble_function(Moose::assemble_matrix);
644 : #else
645 : mooseAssert(
646 : numNonlinearSystems() == 1,
647 : "We should have errored during construction if we had more than one nonlinear system");
648 : mooseAssert(numLinearSystems() == 0,
649 : "We should have errored during construction if we had any linear systems");
650 554 : if (isNonlinearEigenvalueSolver(0))
651 : // We don't need to assemble before the solve
652 486 : _nl_eigen->sys().assemble_before_solve = false;
653 : else
654 68 : _nl_eigen->sys().attach_assemble_function(Moose::assemble_matrix);
655 :
656 : // If matrix_free=true, this tells Libmesh to use shell matrices
657 1028 : _nl_eigen->sys().use_shell_matrices(solverParams(0)._eigen_matrix_free &&
658 474 : !solverParams(0)._eigen_matrix_vector_mult);
659 : // We need to tell libMesh if we are using a shell preconditioning matrix
660 554 : _nl_eigen->sys().use_shell_precond_matrix(solverParams(0)._precond_matrix_free);
661 : #endif
662 :
663 554 : FEProblemBase::init();
664 554 : }
665 :
666 : bool
667 704 : EigenProblem::solverSystemConverged(unsigned int)
668 : {
669 704 : if (_solve)
670 704 : return _nl_eigen->converged();
671 : else
672 0 : return true;
673 : }
674 :
675 : bool
676 12195 : EigenProblem::isNonlinearEigenvalueSolver(const unsigned int eigen_sys_num) const
677 : {
678 12195 : const auto & solver_params = solverParams(eigen_sys_num);
679 24300 : return solver_params._eigen_solve_type == Moose::EST_NONLINEAR_POWER ||
680 12105 : solver_params._eigen_solve_type == Moose::EST_NEWTON ||
681 11717 : solver_params._eigen_solve_type == Moose::EST_PJFNK ||
682 30555 : solver_params._eigen_solve_type == Moose::EST_JFNK ||
683 18450 : solver_params._eigen_solve_type == Moose::EST_PJFNKMO;
684 : }
685 :
686 : void
687 1178 : EigenProblem::initPetscOutputAndSomeSolverSettings()
688 : {
689 1178 : _app.getOutputWarehouse().solveSetup();
690 1178 : }
691 :
692 : Real
693 2345 : EigenProblem::formNorm()
694 : {
695 : mooseAssert(_bx_norm_name,
696 : "We should not get here unless a bx_norm postprocessor has been provided");
697 2345 : return getPostprocessorValueByName(*_bx_norm_name);
698 : }
699 : #endif
700 :
701 : std::string
702 542 : EigenProblem::solverTypeString(const unsigned int solver_sys_num)
703 : {
704 542 : return Moose::stringify(solverParams(solver_sys_num)._eigen_solve_type);
705 : }
|
