Merge pull request #27276 from lindsayad/condensed-eigen-system-27262

Allow solution of eigen problems with constraints

framework/include/problems/EigenProblem.h

@@ -128,7 +128,7 @@ class EigenProblem : public FEProblemBase
    * Form several matrices simultaneously
    */
   void computeMatricesTags(const NumericVector<Number> & soln,
-                           const std::vector<std::unique_ptr<SparseMatrix<Number>>> & jacobians,
+                           const std::vector<SparseMatrix<Number> *> & jacobians,
                            const std::set<TagID> & tags);
 
   /**

framework/include/systems/NonlinearEigenSystem.h

@@ -14,7 +14,7 @@
 #include "NonlinearSystemBase.h"
 #include "SlepcEigenSolverConfiguration.h"
 
-#include "libmesh/eigen_system.h"
+#include "libmesh/condensed_eigen_system.h"
 
 // forward declarations
 class EigenProblem;
@@ -23,6 +23,11 @@ class ResidualObject;
 
 #ifdef LIBMESH_HAVE_SLEPC
 
+namespace Moose
+{
+void assemble_matrix(EquationSystems & es, const std::string & system_name);
+}
+
 /**
  * Nonlinear eigenvalue system to be solved
  */
@@ -87,7 +92,7 @@ class NonlinearEigenSystem : public NonlinearSystemBase
    */
   virtual EPS getEPS();
 
-  EigenSystem & sys() { return _eigen_sys; }
+  CondensedEigenSystem & sys() { return _eigen_sys; }
 
   /**
    * For eigenvalue problems (including standard and generalized), inhomogeneous (Dirichlet or
@@ -168,13 +173,22 @@ 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;
+
 protected:
   virtual void postAddResidualObject(ResidualObject & object) override;
 
   void computeScalingJacobian() override;
   void computeScalingResidual() override;
 
-  EigenSystem & _eigen_sys;
+  CondensedEigenSystem & _eigen_sys;
   EigenProblem & _eigen_problem;
   std::unique_ptr<SlepcEigenSolverConfiguration> _solver_configuration;
   std::vector<std::pair<Real, Real>> _eigen_values;
@@ -189,6 +203,10 @@ class NonlinearEigenSystem : public NonlinearSystemBase
   bool _precond_matrix_includes_eigen;
   // Libmesh preconditioner
   Preconditioner<Number> * _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

framework/include/systems/NonlinearSystem.h

@@ -32,7 +32,7 @@ class NonlinearSystem : public NonlinearSystemBase, public NonlinearImplicitSyst
 
   virtual void solve() override;
 
-  void init() override;
+  virtual void preInit() override;
 
   /**
    * Quit the current solve as soon as possible.

framework/include/systems/NonlinearSystemBase.h

@@ -70,7 +70,7 @@ class NonlinearSystemBase : public SolverSystem, public PerfGraphInterface
   NonlinearSystemBase(FEProblemBase & problem, System & sys, const std::string & name);
   virtual ~NonlinearSystemBase();
 
-  virtual void init() override;
+  virtual void preInit() override;
 
   bool computedScalingJacobian() const { return _computed_scaling; }
 

framework/include/systems/SolverSystem.h

@@ -28,7 +28,7 @@ class SolverSystem : public SystemBase
                Moose::VarKindType var_kind);
   virtual ~SolverSystem();
 
-  virtual void init() override;
+  virtual void preInit() override;
   virtual void restoreSolutions() override final;
 
   void serializeSolution();

framework/include/systems/SystemBase.h

@@ -151,19 +151,32 @@ class SystemBase : public libMesh::ParallelObject, public ConsoleStreamInterface
   virtual const System & system() const = 0;
 
   /**
-   * Initialize the system
+   * This is called prior to the libMesh system has been init'd. MOOSE system wrappers can use this
+   * method to add vectors and matrices to the libMesh system
    */
-  virtual void init(){};
+  virtual void preInit() {}
+
+  /*
+   * This is called after the libMesh system has been init'd. This can be used to initialize MOOSE
+   * system data that relies on the libMesh system data already being initialized
+   */
+  virtual void postInit() {}
+
+  /**
+   * Reinitialize the system when the degrees of freedom in this system have changed. This is called
+   * after the libMesh system has been reinit'd
+   */
+  virtual void reinit() {}
 
   /**
    * Called only once, just before the solve begins so objects can do some precalculations
    */
-  virtual void initializeObjects(){};
+  virtual void initializeObjects() {}
 
   /**
    * Update the system (doing libMesh magic)
    */
-  virtual void update(bool update_libmesh_system = true);
+  void update();
 
   /**
    * Solve the system (using libMesh magic)

framework/include/utils/SlepcSupport.h

@@ -79,14 +79,32 @@ void clearFreeNonlinearPowerIterations(const InputParameters & params);
 
 /**
  * Form matrix according to tag
- */
-void moosePetscSNESFormMatrixTag(SNES snes, Vec x, Mat mat, void * ctx, TagID tag);
-
-/**
- * Form multiple matrices for multiple tags
- */
-void moosePetscSNESFormMatricesTags(
-    SNES snes, Vec x, std::vector<Mat> & mats, void * ctx, const std::set<TagID> & tags);
+ * @param snes The PETSc nonlinear solver context associated with the eigensolver
+ * @param x The current eigen solution vector from SLEPc. This vector will have a size corresponding
+ * to the number of unconstrained degrees of freedom
+ * @param eigen_mat The matrix we want to form coming from SLEPc. As for \p x, this will be sized
+ * corresponding to the number of unconstrained degrees of freedom
+ * @param all_dofs_mat This matrix corresponds to the "global" matrix, representing both
+ * unconstrained and constrained degrees of freedom. We will pass this matrix to MOOSE's assembly
+ * framework and then if the problem has constraints we will selectively copy the submatrix
+ * representing unconstrained dofs to \p eigen_mat. (Else this matrix and \p eigen_mat are the same
+ * so no copy is required)
+ * @param ctx The eigensolver context which corresponds to the MOOSE \p EigenProblem
+ * @param tag The MOOSE \p TagID for the matrix
+ */
+void moosePetscSNESFormMatrixTag(
+    SNES snes, Vec x, Mat eigen_mat, SparseMatrix<Number> & all_dofs_mat, void * ctx, TagID tag);
+
+/**
+ * Form multiple matrices for multiple tags. This function is conceptually the same as \p
+ * moosePetscSNESForMatrixTag() but allows passing in multiple matrices to be formed
+ */
+void moosePetscSNESFormMatricesTags(SNES snes,
+                                    Vec x,
+                                    std::vector<Mat> & eigen_mats,
+                                    std::vector<SparseMatrix<Number> *> & all_dofs_mats,
+                                    void * ctx,
+                                    const std::set<TagID> & tags);
 
 /**
  * Form Jacobian matrix A. It is a SLEPc callback

framework/src/problems/DisplacedProblem.C

@@ -164,21 +164,20 @@ DisplacedProblem::init()
   for (auto & nl : _displaced_solver_systems)
   {
     nl->dofMap().attach_extra_send_list_function(&extraSendList, nl.get());
-    nl->init();
+    nl->preInit();
   }
 
   _displaced_aux->dofMap().attach_extra_send_list_function(&extraSendList, _displaced_aux.get());
-  _displaced_aux->init();
+  _displaced_aux->preInit();
 
   {
     TIME_SECTION("eq::init", 2, "Initializing Displaced Equation System");
     _eq.init();
   }
 
-  /// Get face types properly set for variables
   for (auto & nl : _displaced_solver_systems)
-    nl->update(/*update_libmesh_system=*/false);
-  _displaced_aux->update(/*update_libmesh_system=*/false);
+    nl->postInit();
+  _displaced_aux->postInit();
 
   _mesh.meshChanged();
 
@@ -292,10 +291,6 @@ DisplacedProblem::updateMesh(bool mesh_changing)
   if (haveFV())
     _mesh.setupFiniteVolumeMeshData();
 
-  for (auto & disp_nl : _displaced_solver_systems)
-    disp_nl->update(false);
-  _displaced_aux->update(false);
-
   // Update the geometric searches that depend on the displaced mesh. This call can end up running
   // NearestNodeThread::operator() which has a throw inside of it. We need to catch it and make sure
   // it's propagated to all processes before updating the point locator because the latter requires
@@ -1104,6 +1099,11 @@ DisplacedProblem::meshChanged()
   // EquationSystems::reinit only prolongs/restricts the solution vectors, which is something that
   // needs to happen for every step of mesh adaptivity.
   _eq.reinit();
+  // Since the mesh has changed, we need to make sure that we update any of our
+  // MOOSE-system specific data.
+  for (auto & nl : _displaced_solver_systems)
+    nl->reinit();
+  _displaced_aux->reinit();
 
   // We've performed some mesh adaptivity. We need to
   // clear any quadrature nodes such that when we build the boundary node lists in
@@ -1116,12 +1116,6 @@ DisplacedProblem::meshChanged()
   // then reinitialize GeometricSearchData such that we have all the correct geometric information
   // for the changed mesh
   updateMesh(/*mesh_changing=*/true);
-
-  // Since the mesh has changed, we need to make sure that we update any of our
-  // MOOSE-system specific data. libmesh system data has already been updated
-  for (auto & nl : _displaced_solver_systems)
-    nl->update(/*update_libmesh_system=*/false);
-  _displaced_aux->update(/*update_libmesh_system=*/false);
 }
 
 void

framework/src/problems/EigenProblem.C

@@ -194,10 +194,9 @@ EigenProblem::computeJacobianTag(const NumericVector<Number> & soln,
 }
 
 void
-EigenProblem::computeMatricesTags(
-    const NumericVector<Number> & soln,
-    const std::vector<std::unique_ptr<SparseMatrix<Number>>> & jacobians,
-    const std::set<TagID> & tags)
+EigenProblem::computeMatricesTags(const NumericVector<Number> & soln,
+                                  const std::vector<SparseMatrix<Number> *> & jacobians,
+                                  const std::set<TagID> & tags)
 {
   TIME_SECTION("computeMatricesTags", 3);
 
@@ -622,7 +621,17 @@ EigenProblem::setNormalization(const PostprocessorName & pp, const Real value)
 void
 EigenProblem::init()
 {
-#if !PETSC_RELEASE_LESS_THAN(3, 13, 0)
+#if PETSC_RELEASE_LESS_THAN(3, 13, 0)
+  // Prior to Slepc 3.13 we did not have a nonlinear eigenvalue solver so we must always assemble
+  // before the solve
+  _nl_eigen->sys().attach_assemble_function(Moose::assemble_matrix);
+#else
+  if (isNonlinearEigenvalueSolver())
+    // We don't need to assemble before the solve
+    _nl_eigen->sys().assemble_before_solve = false;
+  else
+    _nl_eigen->sys().attach_assemble_function(Moose::assemble_matrix);
+
   // If matrix_free=true, this tells Libmesh to use shell matrices
   _nl_eigen->sys().use_shell_matrices(solverParams()._eigen_matrix_free &&
                                       !solverParams()._eigen_matrix_vector_mult);

framework/src/problems/FEProblemBase.C

@@ -5894,8 +5894,8 @@ FEProblemBase::init()
       nl->turnOffJacobian();
 
   for (auto & sys : _solver_systems)
-    sys->init();
-  _aux->init();
+    sys->preInit();
+  _aux->preInit();
 
   // Build the mortar segment meshes, if they haven't been already, for a couple reasons:
   // 1) Get the ghosting correct for both static and dynamic meshes
@@ -5913,6 +5913,10 @@ FEProblemBase::init()
     es().init();
   }
 
+  for (auto & sys : _solver_systems)
+    sys->postInit();
+  _aux->postInit();
+
   // Now that the equation system and the dof distribution is done, we can generate the
   // finite volume-related parts if needed.
   if (haveFV())
@@ -7670,7 +7674,14 @@ FEProblemBase::meshChangedHelper(bool intermediate_change)
   if (intermediate_change)
     es().reinit_solutions();
   else
+  {
     es().reinit();
+    // Since the mesh has changed, we need to make sure that we update any of our
+    // MOOSE-system specific data.
+    for (auto & sys : _solver_systems)
+      sys->reinit();
+    _aux->reinit();
+  }
 
   // Updating MooseMesh first breaks other adaptivity code, unless we
   // then *again* update the MooseMesh caches.  E.g. the definition of
@@ -7768,12 +7779,6 @@ FEProblemBase::meshChangedHelper(bool intermediate_change)
     setVariableAllDoFMap(_uo_jacobian_moose_vars[0]);
 
   _has_jacobian = false; // we have to recompute jacobian when mesh changed
-
-  // Since the mesh has changed, we need to make sure that we update any of our
-  // MOOSE-system specific data. libmesh system data has already been updated
-  for (auto & sys : _solver_systems)
-    sys->update(/*update_libmesh_system=*/false);
-  _aux->update(/*update_libmesh_system=*/false);
 }
 
 void