diff --git a/framework/include/systems/NonlinearEigenSystem.h b/framework/include/systems/NonlinearEigenSystem.h index 9f235f8432c2..014fa36ae038 100644 --- a/framework/include/systems/NonlinearEigenSystem.h +++ b/framework/include/systems/NonlinearEigenSystem.h @@ -173,6 +173,12 @@ class NonlinearEigenSystem : public NonlinearSystemBase void residualAndJacobianTogether() override; + /** + * Initialize the condensed matrices. This is a no-op if there are no + * constraints in the DofMap + */ + void initializeCondensedMatrices(); + virtual void postInit() override; virtual void reinit() override; @@ -197,6 +203,10 @@ class NonlinearEigenSystem : public NonlinearSystemBase bool _precond_matrix_includes_eigen; // Libmesh preconditioner Preconditioner * _preconditioner; + + /// The number of degrees of freedom constrained at the libMesh level, e.g. via hanging node or + /// periodic boundary constraints + dof_id_type _num_constrained_dofs; }; #else diff --git a/framework/src/systems/NonlinearEigenSystem.C b/framework/src/systems/NonlinearEigenSystem.C index 95c4ee6d2a34..78bd527763f9 100644 --- a/framework/src/systems/NonlinearEigenSystem.C +++ b/framework/src/systems/NonlinearEigenSystem.C @@ -113,7 +113,8 @@ NonlinearEigenSystem::NonlinearEigenSystem(EigenProblem & eigen_problem, const s _work_rhs_vector_AX(addVector("work_rhs_vector_Ax", false, PARALLEL)), _work_rhs_vector_BX(addVector("work_rhs_vector_Bx", false, PARALLEL)), _precond_matrix_includes_eigen(false), - _preconditioner(nullptr) + _preconditioner(nullptr), + _num_constrained_dofs(0) { SlepcEigenSolver * solver = cast_ptr *>(_eigen_sys.eigen_solver.get()); @@ -190,18 +191,46 @@ NonlinearEigenSystem::postAddResidualObject(ResidualObject & object) } } +void +NonlinearEigenSystem::initializeCondensedMatrices() +{ + if (!(_num_constrained_dofs = dofMap().n_constrained_dofs())) + return; + + _eigen_sys.initialize_condensed_dofs(); + const auto m = cast_int(_eigen_sys.local_non_condensed_dofs_vector.size()); + auto M = m; + _communicator.sum(M); + if (_eigen_sys.has_condensed_matrix_A()) + { + _eigen_sys.get_condensed_matrix_A().init(M, M, m, m); + // A bit ludicrously MatCopy requires the matrix being copied to to be assembled + _eigen_sys.get_condensed_matrix_A().close(); + } + if (_eigen_sys.has_condensed_matrix_B()) + { + _eigen_sys.get_condensed_matrix_B().init(M, M, m, m); + _eigen_sys.get_condensed_matrix_B().close(); + } + if (_eigen_sys.has_condensed_precond_matrix()) + { + _eigen_sys.get_condensed_precond_matrix().init(M, M, m, m); + _eigen_sys.get_condensed_precond_matrix().close(); + } +} + void NonlinearEigenSystem::postInit() { NonlinearSystemBase::postInit(); - _eigen_sys.initialize_condensed_matrices(); + initializeCondensedMatrices(); } void NonlinearEigenSystem::reinit() { NonlinearSystemBase::reinit(); - _eigen_sys.initialize_condensed_matrices(); + initializeCondensedMatrices(); } void @@ -216,7 +245,7 @@ NonlinearEigenSystem::solve() // We apply initial guess for only nonlinear solver if (_eigen_problem.isNonlinearEigenvalueSolver()) { - if (_eigen_sys.have_condensed_dofs()) + if (_num_constrained_dofs) { subvec = solution().get_subvector(_eigen_sys.local_non_condensed_dofs_vector); _eigen_sys.set_initial_space(*subvec); @@ -251,7 +280,7 @@ NonlinearEigenSystem::solve() if (n_converged_eigenvalues) getConvergedEigenpair(_eigen_problem.activeEigenvalueIndex()); - if (_eigen_problem.isNonlinearEigenvalueSolver() && _eigen_sys.have_condensed_dofs()) + if (_eigen_problem.isNonlinearEigenvalueSolver() && _num_constrained_dofs) solution().restore_subvector(std::move(subvec), _eigen_sys.local_non_condensed_dofs_vector); } @@ -261,7 +290,7 @@ NonlinearEigenSystem::attachSLEPcCallbacks() // Tell libmesh not to close matrices before solve _eigen_sys.get_eigen_solver().set_close_matrix_before_solve(false); - if (_eigen_sys.have_condensed_dofs()) + if (_num_constrained_dofs) { // Condensed Matrix A if (_eigen_sys.has_condensed_matrix_A())