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 4191 : EigenProblem::validParams()
36 : {
37 4191 : InputParameters params = FEProblemBase::validParams();
38 8382 : params.addClassDescription("Problem object for solving an eigenvalue problem.");
39 12573 : params.addParam<bool>("negative_sign_eigen_kernel",
40 8382 : 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 12573 : params.addParam<unsigned int>(
46 : "active_eigen_index",
47 8382 : 0,
48 : "Which eigenvector is used to compute residual and also associated to nonlinear variable");
49 16764 : params.addParam<PostprocessorName>("bx_norm", "A postprocessor describing the norm of Bx");
50 :
51 12573 : params.addParamNamesToGroup("negative_sign_eigen_kernel active_eigen_index bx_norm",
52 : "Eigenvalue solve");
53 :
54 4191 : return params;
55 0 : }
56 :
57 565 : 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 565 : _n_eigen_pairs_required(1),
63 565 : _generalized_eigenvalue_problem(false),
64 1130 : _negative_sign_eigen_kernel(getParam<bool>("negative_sign_eigen_kernel")),
65 1130 : _active_eigen_index(getParam<unsigned int>("active_eigen_index")),
66 565 : _do_free_power_iteration(false),
67 565 : _output_inverse_eigenvalue(false),
68 565 : _on_linear_solver(false),
69 565 : _matrices_formed(false),
70 565 : _constant_matrices(false),
71 565 : _has_normalization(false),
72 565 : _normal_factor(1.0),
73 1130 : _first_solve(declareRestartableData<bool>("first_solve", true)),
74 1761 : _bx_norm_name(isParamValid("bx_norm")
75 565 : ? std::make_optional(getParam<PostprocessorName>("bx_norm"))
76 1130 : : std::nullopt)
77 : #endif
78 : {
79 : #ifdef LIBMESH_HAVE_SLEPC
80 565 : if (_nl_sys_names.size() > 1)
81 0 : paramError("nl_sys_names",
82 : "eigen problems do not currently support multiple nonlinear eigen systems");
83 565 : if (_linear_sys_names.size())
84 0 : paramError("linear_sys_names", "EigenProblem only works with a single nonlinear eigen system");
85 :
86 1130 : for (const auto i : index_range(_nl_sys_names))
87 : {
88 565 : const auto & sys_name = _nl_sys_names[i];
89 565 : auto & nl = _nl[i];
90 565 : nl = std::make_shared<NonlinearEigenSystem>(*this, sys_name);
91 565 : _nl_eigen = std::dynamic_pointer_cast<NonlinearEigenSystem>(nl);
92 565 : _current_nl_sys = nl.get();
93 565 : _solver_systems[i] = std::dynamic_pointer_cast<SolverSystem>(nl);
94 565 : nl->system().prefer_hash_table_matrix_assembly(_use_hash_table_matrix_assembly);
95 : }
96 :
97 565 : _aux = std::make_shared<AuxiliarySystem>(*this, "aux0");
98 :
99 565 : newAssemblyArray(_solver_systems);
100 :
101 565 : FEProblemBase::initNullSpaceVectors(parameters, _nl);
102 :
103 1130 : es().parameters.set<EigenProblem *>("_eigen_problem") = this;
104 : #else
105 : mooseError("Need to install SLEPc to solve eigenvalue problems, please reconfigure libMesh\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 if any
114 565 : createTagVectors();
115 :
116 : // Create extra solution vectors if any
117 565 : createTagSolutions();
118 565 : }
119 :
120 : #ifdef LIBMESH_HAVE_SLEPC
121 : void
122 562 : EigenProblem::setEigenproblemType(Moose::EigenProblemType eigen_problem_type)
123 : {
124 562 : switch (eigen_problem_type)
125 : {
126 27 : case Moose::EPT_HERMITIAN:
127 27 : _nl_eigen->sys().set_eigenproblem_type(libMesh::HEP);
128 27 : _generalized_eigenvalue_problem = false;
129 27 : 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 523 : case Moose::EPT_GEN_NON_HERMITIAN:
147 523 : _nl_eigen->sys().set_eigenproblem_type(libMesh::GNHEP);
148 523 : _generalized_eigenvalue_problem = true;
149 523 : 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 562 : }
163 :
164 : void
165 4827 : EigenProblem::execute(const ExecFlagType & exec_type)
166 : {
167 4827 : 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 511 : preScaleEigenVector(std::pair<Real, Real>(_initial_eigenvalue, 0));
171 :
172 4827 : FEProblemBase::execute(exec_type);
173 4827 : }
174 :
175 : void
176 3567 : EigenProblem::computeJacobianTag(const NumericVector<Number> & soln,
177 : SparseMatrix<Number> & jacobian,
178 : TagID tag)
179 : {
180 10701 : 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 3567 : _nl_eigen->disassociateDefaultMatrixTags();
185 :
186 : // Clear FE tags and first add the specific tag associated with the Jacobian
187 3567 : _fe_matrix_tags.clear();
188 3567 : _fe_matrix_tags.insert(tag);
189 :
190 : // Add any other user-added matrix tags if they have associated matrices
191 3567 : const auto & matrix_tags = getMatrixTags();
192 21844 : for (const auto & matrix_tag : matrix_tags)
193 18277 : if (_nl_eigen->hasMatrix(matrix_tag.second))
194 442 : _fe_matrix_tags.insert(matrix_tag.second);
195 :
196 3567 : _nl_eigen->setSolution(soln);
197 :
198 3567 : _nl_eigen->associateMatrixToTag(jacobian, tag);
199 :
200 3567 : setCurrentNonlinearSystem(_nl_eigen->number());
201 3567 : computeJacobianTags(_fe_matrix_tags);
202 :
203 3567 : _nl_eigen->disassociateMatrixFromTag(jacobian, tag);
204 3567 : }
205 :
206 : void
207 348 : EigenProblem::computeMatricesTags(const NumericVector<Number> & soln,
208 : const std::vector<SparseMatrix<Number> *> & jacobians,
209 : const std::set<TagID> & tags)
210 : {
211 1044 : TIME_SECTION("computeMatricesTags", 3);
212 :
213 348 : 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 348 : _nl_eigen->disassociateDefaultMatrixTags();
222 :
223 348 : _fe_matrix_tags.clear();
224 :
225 348 : _nl_eigen->setSolution(soln);
226 :
227 348 : unsigned int i = 0;
228 1044 : for (auto tag : tags)
229 696 : _nl_eigen->associateMatrixToTag(*(jacobians[i++]), tag);
230 :
231 348 : setCurrentNonlinearSystem(_nl_eigen->number());
232 348 : computeJacobianTags(tags);
233 :
234 348 : i = 0;
235 1044 : for (auto tag : tags)
236 696 : _nl_eigen->disassociateMatrixFromTag(*(jacobians[i++]), tag);
237 348 : }
238 :
239 : void
240 324 : EigenProblem::computeJacobianBlocks(std::vector<JacobianBlock *> & blocks,
241 : const unsigned int nl_sys_num)
242 : {
243 972 : TIME_SECTION("computeJacobianBlocks", 3);
244 324 : setCurrentNonlinearSystem(nl_sys_num);
245 :
246 324 : if (_displaced_problem)
247 0 : computeSystems(EXEC_PRE_DISPLACE);
248 :
249 324 : computeSystems(EXEC_NONLINEAR);
250 :
251 324 : _currently_computing_jacobian = true;
252 :
253 648 : _current_nl_sys->computeJacobianBlocks(blocks, {_nl_eigen->precondMatrixTag()});
254 :
255 324 : _currently_computing_jacobian = false;
256 324 : }
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 66 : 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 17800 : EigenProblem::computeResidualTag(const NumericVector<Number> & soln,
296 : NumericVector<Number> & residual,
297 : TagID tag)
298 : {
299 53400 : 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 17800 : _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 17800 : _fe_vector_tags.insert(tag);
308 :
309 : // Add any other user-added vector residual tags if they have associated vectors
310 17800 : const auto & residual_vector_tags = getVectorTags(Moose::VECTOR_TAG_RESIDUAL);
311 89330 : for (const auto & vector_tag : residual_vector_tags)
312 71530 : if (_nl_eigen->hasVector(vector_tag._id))
313 330 : _fe_vector_tags.insert(vector_tag._id);
314 :
315 17800 : _nl_eigen->associateVectorToTag(residual, tag);
316 :
317 17800 : _nl_eigen->setSolution(soln);
318 :
319 17800 : setCurrentNonlinearSystem(_nl_eigen->number());
320 17800 : computeResidualTags(_fe_vector_tags);
321 17800 : _fe_vector_tags.clear();
322 :
323 17800 : _nl_eigen->disassociateVectorFromTag(residual, tag);
324 17800 : }
325 :
326 : void
327 12815 : EigenProblem::computeResidualAB(const NumericVector<Number> & soln,
328 : NumericVector<Number> & residualA,
329 : NumericVector<Number> & residualB,
330 : TagID tagA,
331 : TagID tagB)
332 : {
333 38445 : 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 12815 : _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 12815 : _fe_vector_tags.insert(tagA);
342 12815 : _fe_vector_tags.insert(tagB);
343 :
344 : // Add any other user-added vector residual tags if they have associated vectors
345 12815 : const auto & residual_vector_tags = getVectorTags(Moose::VECTOR_TAG_RESIDUAL);
346 64451 : for (const auto & vector_tag : residual_vector_tags)
347 51636 : if (_nl_eigen->hasVector(vector_tag._id))
348 376 : _fe_vector_tags.insert(vector_tag._id);
349 :
350 12815 : _nl_eigen->associateVectorToTag(residualA, tagA);
351 12815 : _nl_eigen->associateVectorToTag(residualB, tagB);
352 :
353 12815 : _nl_eigen->setSolution(soln);
354 :
355 12815 : computeResidualTags(_fe_vector_tags);
356 12815 : _fe_vector_tags.clear();
357 :
358 12815 : _nl_eigen->disassociateVectorFromTag(residualA, tagA);
359 12815 : _nl_eigen->disassociateVectorFromTag(residualB, tagB);
360 12815 : }
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 1061 : EigenProblem::adjustEigenVector(const Real value, bool scaling)
390 : {
391 1061 : std::vector<VariableName> var_names = getVariableNames();
392 3105 : for (auto & vn : var_names)
393 : {
394 2044 : MooseVariableBase * var = nullptr;
395 2044 : if (hasScalarVariable(vn))
396 46 : var = &getScalarVariable(0, vn);
397 : else
398 1998 : var = &getVariable(0, vn);
399 : // Do operations for only eigen variable
400 2044 : if (var->eigen())
401 2656 : for (unsigned int vc = 0; vc < var->count(); ++vc)
402 : {
403 1433 : std::set<dof_id_type> var_indices;
404 1433 : _nl_eigen->system().local_dof_indices(var->number() + vc, var_indices);
405 90910 : for (const auto & dof : var_indices)
406 89477 : _nl_eigen->solution().set(dof, scaling ? (_nl_eigen->solution()(dof) * value) : value);
407 1433 : }
408 : }
409 :
410 1061 : _nl_eigen->solution().close();
411 1061 : _nl_eigen->update();
412 1061 : }
413 :
414 : void
415 520 : EigenProblem::scaleEigenvector(const Real scaling_factor)
416 : {
417 520 : adjustEigenVector(scaling_factor, true);
418 520 : }
419 :
420 : void
421 541 : 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 541 : adjustEigenVector(initial_value, false);
430 541 : }
431 :
432 : void
433 709 : 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 1418 : computeResidualTag(
437 2127 : *_nl_eigen->currentSolution(), _nl_eigen->residualVectorBX(), _nl_eigen->eigenVectorTag());
438 :
439 : // Eigenvalue magnitude
440 709 : Real v = std::sqrt(eig.first * eig.first + eig.second * eig.second);
441 : // Scaling factor
442 709 : Real factor = 1 / v / (bxNormProvided() ? formNorm() : _nl_eigen->residualVectorBX().l2_norm());
443 : // Scale eigenvector
444 709 : if (!MooseUtils::absoluteFuzzyEqual(factor, 1))
445 511 : scaleEigenvector(factor);
446 709 : }
447 :
448 : void
449 709 : EigenProblem::postScaleEigenVector()
450 : {
451 709 : if (_has_normalization)
452 : {
453 : Real v;
454 9 : 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 9 : v = _normal_factor;
469 :
470 9 : 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 9 : unsigned int itr = 0;
477 :
478 18 : while (!MooseUtils::relativeFuzzyEqual(v, c))
479 : {
480 : // If postprocessor is not defined on eigen variables, scaling might not work
481 9 : 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 9 : scaleEigenvector(v / c);
489 :
490 : // update all aux variables and user objects on linear
491 9 : execute(EXEC_LINEAR);
492 :
493 9 : c = getPostprocessorValueByName(_normalization);
494 :
495 9 : itr++;
496 : }
497 : }
498 709 : }
499 :
500 : void
501 559 : EigenProblem::checkProblemIntegrity()
502 : {
503 559 : FEProblemBase::checkProblemIntegrity();
504 559 : _nl_eigen->checkIntegrity();
505 553 : if (_bx_norm_name)
506 : {
507 33 : if (!isNonlinearEigenvalueSolver(0))
508 0 : paramWarning("bx_norm", "This parameter is only used for nonlinear solve types");
509 33 : else if (auto & pp = getUserObjectBase(_bx_norm_name.value());
510 33 : !pp.getExecuteOnEnum().contains(EXEC_LINEAR))
511 6 : pp.paramError("execute_on",
512 : "If providing the Bx norm, this postprocessor must execute on linear e.g. "
513 : "during residual evaluations");
514 : }
515 550 : }
516 :
517 : void
518 535 : EigenProblem::doFreeNonlinearPowerIterations(unsigned int free_power_iterations)
519 : {
520 : mooseAssert(_current_nl_sys, "This needs to be non-null");
521 :
522 535 : doFreePowerIteration(true);
523 : // Set free power iterations
524 535 : Moose::SlepcSupport::setFreeNonlinearPowerIterations(free_power_iterations);
525 :
526 : // Call solver
527 535 : _current_nl_sys->solve();
528 535 : _current_nl_sys->update();
529 :
530 : // Clear free power iterations
531 535 : auto executioner = getMooseApp().getExecutioner();
532 535 : if (executioner)
533 535 : Moose::SlepcSupport::clearFreeNonlinearPowerIterations(executioner->parameters());
534 : else
535 0 : mooseError("There is no executioner for this moose app");
536 :
537 535 : doFreePowerIteration(false);
538 535 : }
539 :
540 : void
541 718 : 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 718 : if (!_app.isUltimateMaster())
546 132 : LibmeshPetscCall(PetscOptionsPush(_petsc_option_data_base));
547 : #endif
548 :
549 718 : setCurrentNonlinearSystem(nl_sys_num);
550 :
551 718 : if (_solve)
552 : {
553 2127 : 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 709 : _nl_eigen->attachSLEPcCallbacks();
559 :
560 : // If there is an eigenvalue, we scale 1/|Bx| to eigenvalue
561 709 : 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 1363 : if (isNonlinearEigenvalueSolver(nl_sys_num) &&
568 654 : 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 643 : if (solverParams(nl_sys_num)._free_power_iterations && _first_solve)
574 : {
575 517 : _console << std::endl << " -------------------------------" << std::endl;
576 517 : _console << " Free power iteration starts ..." << std::endl;
577 517 : _console << " -------------------------------" << std::endl << std::endl;
578 517 : doFreeNonlinearPowerIterations(solverParams(nl_sys_num)._free_power_iterations);
579 517 : _first_solve = false;
580 : }
581 :
582 : // Let us do extra power iterations here if necessary
583 643 : if (solverParams(nl_sys_num)._extra_power_iterations)
584 : {
585 18 : _console << std::endl << " --------------------------------------" << std::endl;
586 18 : _console << " Extra Free power iteration starts ..." << std::endl;
587 18 : _console << " --------------------------------------" << std::endl << std::endl;
588 18 : doFreeNonlinearPowerIterations(solverParams(nl_sys_num)._extra_power_iterations);
589 : }
590 : }
591 :
592 : // We print this for only nonlinear solver
593 709 : if (isNonlinearEigenvalueSolver(nl_sys_num))
594 : {
595 654 : _console << std::endl << " -------------------------------------" << std::endl;
596 :
597 654 : if (solverParams(nl_sys_num)._eigen_solve_type != Moose::EST_NONLINEAR_POWER)
598 643 : _console << " Nonlinear Newton iteration starts ..." << std::endl;
599 : else
600 11 : _console << " Nonlinear power iteration starts ..." << std::endl;
601 :
602 654 : _console << " -------------------------------------" << std::endl << std::endl;
603 : }
604 :
605 709 : _current_nl_sys->solve();
606 709 : _current_nl_sys->update();
607 :
608 : // with PJFNKMO solve type, we need to evaluate the linear objects to bring them up-to-date
609 709 : if (solverParams(nl_sys_num)._eigen_solve_type == Moose::EST_PJFNKMO)
610 84 : execute(EXEC_LINEAR);
611 :
612 : // Scale eigen vector if users ask
613 709 : postScaleEigenVector();
614 709 : }
615 :
616 : #if !PETSC_RELEASE_LESS_THAN(3, 12, 0)
617 718 : if (!_app.isUltimateMaster())
618 132 : LibmeshPetscCall(PetscOptionsPop());
619 : #endif
620 :
621 : // sync solutions in displaced problem
622 718 : if (_displaced_problem)
623 33 : _displaced_problem->syncSolutions();
624 :
625 : // Reset the matrix flag, so that we reform matrix in next picard iteration
626 718 : _matrices_formed = false;
627 718 : }
628 :
629 : void
630 9 : EigenProblem::setNormalization(const PostprocessorName & pp, const Real value)
631 : {
632 9 : _has_normalization = true;
633 9 : _normalization = pp;
634 9 : _normal_factor = value;
635 9 : }
636 :
637 : void
638 559 : 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 559 : if (isNonlinearEigenvalueSolver(0))
651 : // We don't need to assemble before the solve
652 493 : _nl_eigen->sys().assemble_before_solve = false;
653 : else
654 66 : _nl_eigen->sys().attach_assemble_function(Moose::assemble_matrix);
655 :
656 : // If matrix_free=true, this tells Libmesh to use shell matrices
657 1040 : _nl_eigen->sys().use_shell_matrices(solverParams(0)._eigen_matrix_free &&
658 481 : !solverParams(0)._eigen_matrix_vector_mult);
659 : // We need to tell libMesh if we are using a shell preconditioning matrix
660 559 : _nl_eigen->sys().use_shell_precond_matrix(solverParams(0)._precond_matrix_free);
661 : #endif
662 :
663 559 : FEProblemBase::init();
664 559 : }
665 :
666 : bool
667 709 : EigenProblem::solverSystemConverged(unsigned int)
668 : {
669 709 : if (_solve)
670 709 : return _nl_eigen->converged();
671 : else
672 0 : return true;
673 : }
674 :
675 : bool
676 14649 : EigenProblem::isNonlinearEigenvalueSolver(const unsigned int eigen_sys_num) const
677 : {
678 14649 : const auto & solver_params = solverParams(eigen_sys_num);
679 29186 : return solver_params._eigen_solve_type == Moose::EST_NONLINEAR_POWER ||
680 14537 : solver_params._eigen_solve_type == Moose::EST_NEWTON ||
681 14037 : solver_params._eigen_solve_type == Moose::EST_PJFNK ||
682 35724 : solver_params._eigen_solve_type == Moose::EST_JFNK ||
683 21187 : solver_params._eigen_solve_type == Moose::EST_PJFNKMO;
684 : }
685 :
686 : void
687 1189 : EigenProblem::initPetscOutputAndSomeSolverSettings()
688 : {
689 1189 : _app.getOutputWarehouse().solveSetup();
690 1189 : }
691 :
692 : Real
693 2152 : EigenProblem::formNorm()
694 : {
695 : mooseAssert(_bx_norm_name,
696 : "We should not get here unless a bx_norm postprocessor has been provided");
697 2152 : return getPostprocessorValueByName(*_bx_norm_name);
698 : }
699 : #endif
700 :
701 : std::string
702 550 : EigenProblem::solverTypeString(const unsigned int solver_sys_num)
703 : {
704 550 : return Moose::stringify(solverParams(solver_sys_num)._eigen_solve_type);
705 : }
|
