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 15441 : EigenProblem::validParams()
36 : {
37 15441 : InputParameters params = FEProblemBase::validParams();
38 15441 : params.addClassDescription("Problem object for solving an eigenvalue problem.");
39 46323 : params.addParam<bool>("negative_sign_eigen_kernel",
40 30882 : 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 46323 : params.addParam<unsigned int>(
46 : "active_eigen_index",
47 30882 : 0,
48 : "Which eigenvector is used to compute residual and also associated to nonlinear variable");
49 15441 : params.addParam<PostprocessorName>("bx_norm", "A postprocessor describing the norm of Bx");
50 :
51 15441 : params.addParamNamesToGroup("negative_sign_eigen_kernel active_eigen_index bx_norm",
52 : "Eigenvalue solve");
53 :
54 15441 : return params;
55 0 : }
56 :
57 588 : 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 588 : _n_eigen_pairs_required(1),
63 588 : _generalized_eigenvalue_problem(false),
64 588 : _negative_sign_eigen_kernel(getParam<bool>("negative_sign_eigen_kernel")),
65 588 : _active_eigen_index(getParam<unsigned int>("active_eigen_index")),
66 588 : _do_free_power_iteration(false),
67 588 : _output_inverse_eigenvalue(false),
68 588 : _on_linear_solver(false),
69 588 : _matrices_formed(false),
70 588 : _constant_matrices(false),
71 588 : _has_normalization(false),
72 588 : _normal_factor(1.0),
73 588 : _first_solve(declareRestartableData<bool>("first_solve", true)),
74 1176 : _bx_norm_name(isParamValid("bx_norm")
75 588 : ? std::make_optional(getParam<PostprocessorName>("bx_norm"))
76 1176 : : std::nullopt)
77 : #endif
78 : {
79 : #ifdef LIBMESH_HAVE_SLEPC
80 588 : if (_nl_sys_names.size() > 1)
81 0 : paramError("nl_sys_names",
82 : "eigen problems do not currently support multiple nonlinear eigen systems");
83 588 : if (_linear_sys_names.size())
84 0 : paramError("linear_sys_names", "EigenProblem only works with a single nonlinear eigen system");
85 :
86 1176 : for (const auto i : index_range(_nl_sys_names))
87 : {
88 588 : const auto & sys_name = _nl_sys_names[i];
89 588 : auto & nl = _nl[i];
90 588 : nl = std::make_shared<NonlinearEigenSystem>(*this, sys_name);
91 588 : _nl_eigen = std::dynamic_pointer_cast<NonlinearEigenSystem>(nl);
92 588 : _current_nl_sys = nl.get();
93 588 : _solver_systems[i] = std::dynamic_pointer_cast<SolverSystem>(nl);
94 588 : nl->system().prefer_hash_table_matrix_assembly(_use_hash_table_matrix_assembly);
95 : }
96 :
97 588 : _aux = std::make_shared<AuxiliarySystem>(*this, "aux0");
98 :
99 588 : newAssemblyArray(_solver_systems);
100 :
101 588 : FEProblemBase::initNullSpaceVectors(parameters, _nl);
102 :
103 588 : 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 588 : createTagVectors();
115 :
116 : // Create extra solution vectors if any
117 588 : createTagSolutions();
118 588 : }
119 :
120 : #ifdef LIBMESH_HAVE_SLEPC
121 : void
122 584 : EigenProblem::setEigenproblemType(Moose::EigenProblemType eigen_problem_type)
123 : {
124 584 : switch (eigen_problem_type)
125 : {
126 30 : case Moose::EPT_HERMITIAN:
127 30 : _nl_eigen->sys().set_eigenproblem_type(libMesh::HEP);
128 30 : _generalized_eigenvalue_problem = false;
129 30 : break;
130 :
131 13 : case Moose::EPT_NON_HERMITIAN:
132 13 : _nl_eigen->sys().set_eigenproblem_type(libMesh::NHEP);
133 13 : _generalized_eigenvalue_problem = false;
134 13 : 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 541 : case Moose::EPT_GEN_NON_HERMITIAN:
147 541 : _nl_eigen->sys().set_eigenproblem_type(libMesh::GNHEP);
148 541 : _generalized_eigenvalue_problem = true;
149 541 : 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 584 : }
163 :
164 : void
165 5078 : EigenProblem::execute(const ExecFlagType & exec_type)
166 : {
167 5078 : 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 532 : preScaleEigenVector(std::pair<Real, Real>(_initial_eigenvalue, 0));
171 :
172 5078 : FEProblemBase::execute(exec_type);
173 5078 : }
174 :
175 : void
176 3696 : EigenProblem::computeJacobianTag(const NumericVector<Number> & soln,
177 : SparseMatrix<Number> & jacobian,
178 : TagID tag)
179 : {
180 3696 : 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 3696 : _nl_eigen->disassociateDefaultMatrixTags();
185 :
186 : // Clear FE tags and first add the specific tag associated with the Jacobian
187 3696 : _fe_matrix_tags.clear();
188 3696 : _fe_matrix_tags.insert(tag);
189 :
190 : // Add any other user-added matrix tags if they have associated matrices
191 3696 : const auto & matrix_tags = getMatrixTags();
192 22656 : for (const auto & matrix_tag : matrix_tags)
193 18960 : if (_nl_eigen->hasMatrix(matrix_tag.second))
194 480 : _fe_matrix_tags.insert(matrix_tag.second);
195 :
196 3696 : _nl_eigen->setSolution(soln);
197 :
198 3696 : _nl_eigen->associateMatrixToTag(jacobian, tag);
199 :
200 3696 : setCurrentNonlinearSystem(_nl_eigen->number());
201 3696 : computeJacobianTags(_fe_matrix_tags);
202 :
203 3696 : _nl_eigen->disassociateMatrixFromTag(jacobian, tag);
204 3696 : }
205 :
206 : void
207 380 : EigenProblem::computeMatricesTags(const NumericVector<Number> & soln,
208 : const std::vector<SparseMatrix<Number> *> & jacobians,
209 : const std::set<TagID> & tags)
210 : {
211 380 : TIME_SECTION("computeMatricesTags", 3);
212 :
213 380 : 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 380 : _nl_eigen->disassociateDefaultMatrixTags();
222 :
223 380 : _fe_matrix_tags.clear();
224 :
225 380 : _nl_eigen->setSolution(soln);
226 :
227 380 : unsigned int i = 0;
228 1140 : for (auto tag : tags)
229 760 : _nl_eigen->associateMatrixToTag(*(jacobians[i++]), tag);
230 :
231 380 : setCurrentNonlinearSystem(_nl_eigen->number());
232 380 : computeJacobianTags(tags);
233 :
234 380 : i = 0;
235 1140 : for (auto tag : tags)
236 760 : _nl_eigen->disassociateMatrixFromTag(*(jacobians[i++]), tag);
237 380 : }
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 24 : EigenProblem::computeJacobianAB(const NumericVector<Number> & soln,
260 : SparseMatrix<Number> & jacobianA,
261 : SparseMatrix<Number> & jacobianB,
262 : TagID tagA,
263 : TagID tagB)
264 : {
265 24 : 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 24 : _nl_eigen->disassociateDefaultMatrixTags();
270 :
271 : // Clear FE tags and first add the specific tags associated with the Jacobian
272 24 : _fe_matrix_tags.clear();
273 24 : _fe_matrix_tags.insert(tagA);
274 24 : _fe_matrix_tags.insert(tagB);
275 :
276 : // Add any other user-added matrix tags if they have associated matrices
277 24 : const auto & matrix_tags = getMatrixTags();
278 144 : for (const auto & matrix_tag : matrix_tags)
279 120 : if (_nl_eigen->hasMatrix(matrix_tag.second))
280 0 : _fe_matrix_tags.insert(matrix_tag.second);
281 :
282 24 : _nl_eigen->setSolution(soln);
283 :
284 24 : _nl_eigen->associateMatrixToTag(jacobianA, tagA);
285 24 : _nl_eigen->associateMatrixToTag(jacobianB, tagB);
286 :
287 24 : setCurrentNonlinearSystem(_nl_eigen->number());
288 24 : computeJacobianTags(_fe_matrix_tags);
289 :
290 24 : _nl_eigen->disassociateMatrixFromTag(jacobianA, tagA);
291 24 : _nl_eigen->disassociateMatrixFromTag(jacobianB, tagB);
292 24 : }
293 :
294 : void
295 18452 : EigenProblem::computeResidualTag(const NumericVector<Number> & soln,
296 : NumericVector<Number> & residual,
297 : TagID tag)
298 : {
299 18452 : 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 18452 : _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 18452 : _fe_vector_tags.insert(tag);
308 :
309 : // Add any other user-added vector residual tags if they have associated vectors
310 18452 : const auto & residual_vector_tags = getVectorTags(Moose::VECTOR_TAG_RESIDUAL);
311 92620 : for (const auto & vector_tag : residual_vector_tags)
312 74168 : if (_nl_eigen->hasVector(vector_tag._id))
313 360 : _fe_vector_tags.insert(vector_tag._id);
314 :
315 18452 : _nl_eigen->associateVectorToTag(residual, tag);
316 :
317 18452 : _nl_eigen->setSolution(soln);
318 :
319 18452 : setCurrentNonlinearSystem(_nl_eigen->number());
320 18452 : computeResidualTags(_fe_vector_tags);
321 18452 : _fe_vector_tags.clear();
322 :
323 18452 : _nl_eigen->disassociateVectorFromTag(residual, tag);
324 18452 : }
325 :
326 : void
327 13394 : EigenProblem::computeResidualAB(const NumericVector<Number> & soln,
328 : NumericVector<Number> & residualA,
329 : NumericVector<Number> & residualB,
330 : TagID tagA,
331 : TagID tagB)
332 : {
333 13394 : 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 13394 : _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 13394 : _fe_vector_tags.insert(tagA);
342 13394 : _fe_vector_tags.insert(tagB);
343 :
344 : // Add any other user-added vector residual tags if they have associated vectors
345 13394 : const auto & residual_vector_tags = getVectorTags(Moose::VECTOR_TAG_RESIDUAL);
346 67378 : for (const auto & vector_tag : residual_vector_tags)
347 53984 : if (_nl_eigen->hasVector(vector_tag._id))
348 408 : _fe_vector_tags.insert(vector_tag._id);
349 :
350 13394 : _nl_eigen->associateVectorToTag(residualA, tagA);
351 13394 : _nl_eigen->associateVectorToTag(residualB, tagB);
352 :
353 13394 : _nl_eigen->setSolution(soln);
354 :
355 13394 : computeResidualTags(_fe_vector_tags);
356 13394 : _fe_vector_tags.clear();
357 :
358 13394 : _nl_eigen->disassociateVectorFromTag(residualA, tagA);
359 13394 : _nl_eigen->disassociateVectorFromTag(residualB, tagB);
360 13394 : }
361 :
362 : Real
363 144 : EigenProblem::computeResidualL2Norm()
364 : {
365 576 : computeResidualAB(*_nl_eigen->currentSolution(),
366 144 : _nl_eigen->residualVectorAX(),
367 144 : _nl_eigen->residualVectorBX(),
368 144 : _nl_eigen->nonEigenVectorTag(),
369 144 : _nl_eigen->eigenVectorTag());
370 :
371 144 : Real eigenvalue = 1.0;
372 :
373 144 : if (_active_eigen_index < _nl_eigen->getNumConvergedEigenvalues())
374 120 : eigenvalue = _nl_eigen->getConvergedEigenvalue(_active_eigen_index).first;
375 :
376 : // Scale BX with eigenvalue
377 144 : _nl_eigen->residualVectorBX() *= eigenvalue;
378 :
379 : // Compute entire residual
380 144 : if (_negative_sign_eigen_kernel)
381 144 : _nl_eigen->residualVectorAX() += _nl_eigen->residualVectorBX();
382 : else
383 0 : _nl_eigen->residualVectorAX() -= _nl_eigen->residualVectorBX();
384 :
385 144 : return _nl_eigen->residualVectorAX().l2_norm();
386 : }
387 :
388 : void
389 1100 : EigenProblem::adjustEigenVector(const Real value, bool scaling)
390 : {
391 1100 : std::vector<VariableName> var_names = getVariableNames();
392 3275 : for (auto & vn : var_names)
393 : {
394 2175 : MooseVariableBase * var = nullptr;
395 2175 : if (hasScalarVariable(vn))
396 50 : var = &getScalarVariable(0, vn);
397 : else
398 2125 : var = &getVariable(0, vn);
399 : // Do operations for only eigen variable
400 2175 : if (var->eigen())
401 2790 : for (unsigned int vc = 0; vc < var->count(); ++vc)
402 : {
403 1510 : std::set<dof_id_type> var_indices;
404 1510 : _nl_eigen->system().local_dof_indices(var->number() + vc, var_indices);
405 98594 : for (const auto & dof : var_indices)
406 97084 : _nl_eigen->solution().set(dof, scaling ? (_nl_eigen->solution()(dof) * value) : value);
407 1510 : }
408 : }
409 :
410 1100 : _nl_eigen->solution().close();
411 1100 : _nl_eigen->update();
412 1100 : }
413 :
414 : void
415 542 : EigenProblem::scaleEigenvector(const Real scaling_factor)
416 : {
417 542 : adjustEigenVector(scaling_factor, true);
418 542 : }
419 :
420 : void
421 558 : 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 558 : adjustEigenVector(initial_value, false);
430 558 : }
431 :
432 : void
433 748 : 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 1496 : computeResidualTag(
437 2244 : *_nl_eigen->currentSolution(), _nl_eigen->residualVectorBX(), _nl_eigen->eigenVectorTag());
438 :
439 : // Eigenvalue magnitude
440 748 : Real v = std::sqrt(eig.first * eig.first + eig.second * eig.second);
441 : // Scaling factor
442 748 : Real factor = 1 / v / (bxNormProvided() ? formNorm() : _nl_eigen->residualVectorBX().l2_norm());
443 : // Scale eigenvector
444 748 : if (!MooseUtils::absoluteFuzzyEqual(factor, 1))
445 532 : scaleEigenvector(factor);
446 748 : }
447 :
448 : void
449 748 : EigenProblem::postScaleEigenVector()
450 : {
451 748 : 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 748 : }
499 :
500 : void
501 580 : EigenProblem::checkProblemIntegrity()
502 : {
503 580 : FEProblemBase::checkProblemIntegrity();
504 580 : _nl_eigen->checkIntegrity();
505 572 : if (_bx_norm_name)
506 : {
507 37 : if (!isNonlinearEigenvalueSolver(0))
508 0 : paramWarning("bx_norm", "This parameter is only used for nonlinear solve types");
509 37 : else if (auto & pp = getUserObjectBase(_bx_norm_name.value());
510 37 : !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 568 : }
516 :
517 : void
518 558 : EigenProblem::doFreeNonlinearPowerIterations(unsigned int free_power_iterations)
519 : {
520 : mooseAssert(_current_nl_sys, "This needs to be non-null");
521 :
522 558 : doFreePowerIteration(true);
523 : // Set free power iterations
524 558 : Moose::SlepcSupport::setFreeNonlinearPowerIterations(free_power_iterations);
525 :
526 : // Call solver
527 558 : _current_nl_sys->solve();
528 558 : _current_nl_sys->update();
529 :
530 : // Clear free power iterations
531 558 : auto executioner = getMooseApp().getExecutioner();
532 558 : if (executioner)
533 558 : Moose::SlepcSupport::clearFreeNonlinearPowerIterations(executioner->parameters());
534 : else
535 0 : mooseError("There is no executioner for this moose app");
536 :
537 558 : doFreePowerIteration(false);
538 558 : }
539 :
540 : void
541 758 : 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 758 : if (!_app.isUltimateMaster())
546 144 : LibmeshPetscCall(PetscOptionsPush(_petsc_option_data_base));
547 : #endif
548 :
549 758 : setCurrentNonlinearSystem(nl_sys_num);
550 :
551 758 : if (_solve)
552 : {
553 748 : 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 748 : _nl_eigen->attachSLEPcCallbacks();
559 :
560 : // If there is an eigenvalue, we scale 1/|Bx| to eigenvalue
561 748 : if (_active_eigen_index < _nl_eigen->getNumConvergedEigenvalues())
562 : {
563 216 : std::pair<Real, Real> eig = _nl_eigen->getConvergedEigenvalue(_active_eigen_index);
564 216 : preScaleEigenVector(eig);
565 : }
566 :
567 1436 : if (isNonlinearEigenvalueSolver(nl_sys_num) &&
568 688 : 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 676 : if (solverParams(nl_sys_num)._free_power_iterations && _first_solve)
574 : {
575 538 : _console << std::endl << " -------------------------------" << std::endl;
576 538 : _console << " Free power iteration starts ..." << std::endl;
577 538 : _console << " -------------------------------" << std::endl << std::endl;
578 538 : doFreeNonlinearPowerIterations(solverParams(nl_sys_num)._free_power_iterations);
579 538 : _first_solve = false;
580 : }
581 :
582 : // Let us do extra power iterations here if necessary
583 676 : 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 748 : if (isNonlinearEigenvalueSolver(nl_sys_num))
594 : {
595 688 : _console << std::endl << " -------------------------------------" << std::endl;
596 :
597 688 : if (solverParams(nl_sys_num)._eigen_solve_type != Moose::EST_NONLINEAR_POWER)
598 676 : _console << " Nonlinear Newton iteration starts ..." << std::endl;
599 : else
600 12 : _console << " Nonlinear power iteration starts ..." << std::endl;
601 :
602 688 : _console << " -------------------------------------" << std::endl << std::endl;
603 : }
604 :
605 748 : _current_nl_sys->solve();
606 748 : _current_nl_sys->update();
607 :
608 : // with PJFNKMO solve type, we need to evaluate the linear objects to bring them up-to-date
609 748 : if (solverParams(nl_sys_num)._eigen_solve_type == Moose::EST_PJFNKMO)
610 92 : execute(EXEC_LINEAR);
611 :
612 : // Scale eigen vector if users ask
613 748 : postScaleEigenVector();
614 748 : }
615 :
616 : #if !PETSC_RELEASE_LESS_THAN(3, 12, 0)
617 758 : if (!_app.isUltimateMaster())
618 144 : LibmeshPetscCall(PetscOptionsPop());
619 : #endif
620 :
621 : // sync solutions in displaced problem
622 758 : if (_displaced_problem)
623 36 : _displaced_problem->syncSolutions();
624 :
625 : // Reset the matrix flag, so that we reform matrix in next picard iteration
626 758 : _matrices_formed = false;
627 758 : }
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 580 : 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 580 : if (isNonlinearEigenvalueSolver(0))
651 : // We don't need to assemble before the solve
652 507 : _nl_eigen->sys().assemble_before_solve = false;
653 : else
654 73 : _nl_eigen->sys().attach_assemble_function(Moose::assemble_matrix);
655 :
656 : // If matrix_free=true, this tells Libmesh to use shell matrices
657 1074 : _nl_eigen->sys().use_shell_matrices(solverParams(0)._eigen_matrix_free &&
658 494 : !solverParams(0)._eigen_matrix_vector_mult);
659 : // We need to tell libMesh if we are using a shell preconditioning matrix
660 580 : _nl_eigen->sys().use_shell_precond_matrix(solverParams(0)._precond_matrix_free);
661 : #endif
662 :
663 580 : FEProblemBase::init();
664 580 : }
665 :
666 : bool
667 748 : EigenProblem::solverSystemConverged(unsigned int)
668 : {
669 748 : if (_solve)
670 748 : return _nl_eigen->converged();
671 : else
672 0 : return true;
673 : }
674 :
675 : bool
676 12911 : EigenProblem::isNonlinearEigenvalueSolver(const unsigned int eigen_sys_num) const
677 : {
678 12911 : const auto & solver_params = solverParams(eigen_sys_num);
679 25724 : return solver_params._eigen_solve_type == Moose::EST_NONLINEAR_POWER ||
680 12813 : solver_params._eigen_solve_type == Moose::EST_NEWTON ||
681 12401 : solver_params._eigen_solve_type == Moose::EST_PJFNK ||
682 32375 : solver_params._eigen_solve_type == Moose::EST_JFNK ||
683 19562 : solver_params._eigen_solve_type == Moose::EST_PJFNKMO;
684 : }
685 :
686 : void
687 1246 : EigenProblem::initPetscOutputAndSomeSolverSettings()
688 : {
689 1246 : _app.getOutputWarehouse().solveSetup();
690 1246 : }
691 :
692 : Real
693 2374 : EigenProblem::formNorm()
694 : {
695 : mooseAssert(_bx_norm_name,
696 : "We should not get here unless a bx_norm postprocessor has been provided");
697 2374 : return getPostprocessorValueByName(*_bx_norm_name);
698 : }
699 : #endif
700 :
701 : std::string
702 568 : EigenProblem::solverTypeString(const unsigned int solver_sys_num)
703 : {
704 568 : return Moose::stringify(solverParams(solver_sys_num)._eigen_solve_type);
705 : }
|
