Use tagging system to support eigen kernels in the precond matrix

framework/include/executioners/Eigenvalue.h

@@ -34,6 +34,8 @@ class Eigenvalue : public Steady
 
   Eigenvalue(const InputParameters & parameters);
 
+  virtual void init() override;
+
   virtual void execute() override;
 
   /**

framework/include/problems/EigenProblem.h

@@ -70,11 +70,11 @@ class EigenProblem : public FEProblemBase
                                   TagID tag) override;
 
   /**
-   * Form a Jacobian matrix for all kernels and BCs with given tags
+   * Form several matrices simultaneously
    */
-  virtual void computePrecondMatrixTags(const NumericVector<Number> & soln,
-                                        SparseMatrix<Number> & jacobian,
-                                        std::set<TagID> & tags);
+  void computeMatricesTags(const NumericVector<Number> & soln,
+                           const std::vector<std::unique_ptr<SparseMatrix<Number>>> & jacobians,
+                           const std::set<TagID> & tags);
 
   /**
    * Form two Jacobian matrices, whre each is associateed with one tag, through one

framework/include/systems/NonlinearEigenSystem.h

@@ -70,6 +70,12 @@ class NonlinearEigenSystem : public NonlinearSystemBase
   template <typename T>
   void addEigenTagToMooseObjects(MooseObjectTagWarehouse<T> & warehouse);
 
+  /**
+   * Add the precond tag to eigen kernels
+   */
+  template <typename T>
+  void addPrecondTagToMooseObjects(MooseObjectTagWarehouse<T> & warehouse);
+
   virtual void initialSetup() override;
 
   void attachSLEPcCallbacks();
@@ -168,7 +174,7 @@ class NonlinearEigenSystem : public NonlinearSystemBase
 
   bool precondMatrixIncludesEigenKernels() { return _precond_matrix_includes_eigen; }
 
-  std::set<TagID> & precondMatrixTags() { return _precond_tags; }
+  TagID precondMatrixTag() { return _precond_tag; }
 
   virtual void attachMoosePreconditioner(Preconditioner<Number> * preconditioner) override;
 
@@ -192,9 +198,9 @@ class NonlinearEigenSystem : public NonlinearSystemBase
   TagID _Bx_tag;
   TagID _A_tag;
   TagID _B_tag;
+  TagID _precond_tag;
   bool _precond_matrix_includes_eigen;
   // Tags for computing preconditioning matrix
-  std::set<TagID> _precond_tags;
   Preconditioner<Number> * _moose_preconditioner;
 };
 

framework/include/utils/SlepcSupport.h

@@ -45,9 +45,9 @@ void slepcSetOptions(EigenProblem & eigen_problem, const InputParameters & param
 void setSlepcEigenSolverTolerances(EigenProblem & eigen_problem, const InputParameters & params);
 void setSlepcOutputOptions(EigenProblem & eigen_problem);
 
-void moosePetscSNESFormJacobian(SNES snes, Vec x, Mat jac, Mat pc, void * ctx, TagID tag);
-void moosePetscSNESFormJacobianTags(
-    SNES snes, Vec x, Mat jac, Mat pc, void * ctx, std::set<TagID> & tags);
+void moosePetscSNESFormMatrixTag(SNES snes, Vec x, Mat mat, void * ctx, TagID tag);
+void moosePetscSNESFormMatricesTags(
+    SNES snes, Vec x, std::vector<Mat> & mats, void * ctx, const std::set<TagID> & tags);
 PetscErrorCode mooseSlepcEigenFormJacobianA(SNES snes, Vec x, Mat jac, Mat pc, void * ctx);
 PetscErrorCode mooseSlepcEigenFormJacobianB(SNES snes, Vec x, Mat jac, Mat pc, void * ctx);
 PetscErrorCode mooseSlepcEigenFormFunctionA(SNES snes, Vec x, Vec r, void * ctx);

framework/src/base/Assembly.C

@@ -3890,8 +3890,8 @@ Assembly::setCachedJacobianContributions()
       // First zero the rows (including the diagonals) to prepare for
       // setting the cached values.
       // If we did not cache anything yet, we should skip
-      if (_cached_jacobian_contribution_rows[tag].size())
-        _sys.getMatrix(tag).zero_rows(_cached_jacobian_contribution_rows[tag], 0.0);
+      // if (_cached_jacobian_contribution_rows[tag].size())
+      _sys.getMatrix(tag).zero_rows(_cached_jacobian_contribution_rows[tag], 0.0);
 
       // TODO: Use SparseMatrix::set_values() for efficiency
       for (MooseIndex(_cached_jacobian_contribution_vals) i = 0;

framework/src/executioners/Eigenvalue.C

@@ -80,12 +80,18 @@ Eigenvalue::Eigenvalue(const InputParameters & parameters)
 }
 
 void
-Eigenvalue::execute()
+Eigenvalue::init()
 {
   // Set a flag to nonlinear eigen system
   _eigen_problem.getNonlinearEigenSystem().precondMatrixIncludesEigenKernels(
       getParam<bool>("precond_matrix_includes_eigen"));
 
+  Steady::init();
+}
+
+void
+Eigenvalue::execute()
+{
 #if LIBMESH_HAVE_SLEPC
 #if PETSC_RELEASE_LESS_THAN(3, 12, 0)
   // Make sure the SLEPc options are setup for this app

framework/src/problems/EigenProblem.C

@@ -145,41 +145,34 @@ EigenProblem::computeJacobianTag(const NumericVector<Number> & soln,
 }
 
 void
-EigenProblem::computePrecondMatrixTags(const NumericVector<Number> & soln,
-                                       SparseMatrix<Number> & jacobian,
-                                       std::set<TagID> & tags)
+EigenProblem::computeMatricesTags(
+    const NumericVector<Number> & soln,
+    const std::vector<std::unique_ptr<SparseMatrix<Number>>> & jacobians,
+    const std::set<TagID> & tags)
 {
   TIME_SECTION(_compute_jacobian_tag_timer);
 
+  if (jacobians.size() != tags.size())
+    mooseError("The number of matrices ",
+               jacobians.size(),
+               " does not equal the number of tags ",
+               tags.size());
+
   _fe_matrix_tags.clear();
 
   _nl_eigen->setSolution(soln);
 
   _nl_eigen->disassociateAllTaggedMatrices();
 
+  unsigned int i = 0;
   for (auto tag : tags)
-  {
-    _nl_eigen->associateMatrixToTag(jacobian, tag);
-    // We need to shut off eigen nodal BCs
-    // In "A+B", the corresponding rows of B are just zero
-    // Turning them off avoid overwriting
-    if (tag == _nl_eigen->eigenMatrixTag())
-    {
-      _nl_eigen->turnOffOnEigenNodalBCs(false);
-    }
-  }
+    _nl_eigen->associateMatrixToTag(*(jacobians[i++]), tag);
 
   computeJacobianTags(tags);
 
+  i = 0;
   for (auto tag : tags)
-  {
-    _nl_eigen->disassociateMatrixFromTag(jacobian, tag);
-    // Turn all nodal BCs back
-    if (tag == _nl_eigen->eigenMatrixTag())
-    {
-      _nl_eigen->turnOffOnEigenNodalBCs(true);
-    }
-  }
+    _nl_eigen->disassociateMatrixFromTag(*(jacobians[i++]), tag);
 }
 
 void
@@ -197,7 +190,7 @@ EigenProblem::computeJacobianBlocks(std::vector<JacobianBlock *> & blocks)
 
   _currently_computing_jacobian = true;
 
-  _nl->computeJacobianBlocks(blocks, _nl_eigen->precondMatrixTags());
+  _nl->computeJacobianBlocks(blocks, {_nl_eigen->precondMatrixTag()});
 
   _currently_computing_jacobian = false;
 }

framework/src/systems/NonlinearEigenSystem.C

@@ -51,9 +51,9 @@ assemble_matrix(EquationSystems & es, const std::string & system_name)
   // we only need to form a preconditioning matrix
   if (eigen_system.use_shell_matrices() && !eigen_system.use_shell_precond_matrix())
   {
-    p->computePrecondMatrixTags(*eigen_system.current_local_solution.get(),
-                                *eigen_system.precond_matrix,
-                                eigen_nl.precondMatrixTags());
+    p->computeJacobianTag(*eigen_system.current_local_solution.get(),
+                          *eigen_system.precond_matrix,
+                          eigen_nl.precondMatrixTag());
     return;
   }
 #endif
@@ -101,9 +101,13 @@ NonlinearEigenSystem::NonlinearEigenSystem(EigenProblem & eigen_problem, const s
 
   _A_tag = eigen_problem.addMatrixTag("A_tag");
 
-  _precond_tags.insert(_A_tag);
-
   _B_tag = eigen_problem.addMatrixTag("Eigen");
+
+  // By default, _precond_tag and _A_tag will share the same
+  // objects. If we want to include eigen contributions to
+  // the preconditioning matrix, and then _precond_tag will
+  // point to part of "B" objects
+  _precond_tag = eigen_problem.addMatrixTag("Eigen_precond");
 }
 
 void
@@ -120,6 +124,36 @@ NonlinearEigenSystem::initialSetup()
   addEigenTagToMooseObjects(_scalar_kernels);
   // IntegratedBCs
   addEigenTagToMooseObjects(_integrated_bcs);
+  // If the precond matrix needs to include eigen kernels,
+  // we assign precond tag to all objects except nodal BCs.
+  // Mathematically speaking, we can not add eigen BCs in
+  // since they will overwrite non-eigen nodal BCs contributions.
+  if (_precond_matrix_includes_eigen)
+  {
+    // DG kernels
+    addPrecondTagToMooseObjects(_dg_kernels);
+    // Regular kernels
+    addPrecondTagToMooseObjects(_kernels);
+    // Scalar kernels
+    addPrecondTagToMooseObjects(_scalar_kernels);
+    // IntegratedBCs
+    addPrecondTagToMooseObjects(_integrated_bcs);
+  }
+}
+
+template <typename T>
+void
+NonlinearEigenSystem::addPrecondTagToMooseObjects(MooseObjectTagWarehouse<T> & warehouse)
+{
+  for (THREAD_ID tid = 0; tid < warehouse.numThreads(); tid++)
+  {
+    // Get all objects out from the warehouse
+    auto & objects = warehouse.getObjects(tid);
+
+    // Assign precond tag to all objects
+    for (auto & object : objects)
+      object->useMatrixTag(_precond_tag);
+  }
 }
 
 template <typename T>
@@ -135,13 +169,21 @@ NonlinearEigenSystem::addEigenTagToMooseObjects(MooseObjectTagWarehouse<T> & war
       auto & vtags = object->getVectorTags();
       // If this is not an eigen kernel
       if (vtags.find(_Bx_tag) == vtags.end())
+      {
         object->useVectorTag(_Ax_tag);
+        // Noneigen Kernels
+        object->useMatrixTag(_precond_tag);
+      }
       else // also associate eigen matrix tag if this is a eigen kernel
         object->useMatrixTag(_B_tag);
 
       auto & mtags = object->getMatrixTags();
       if (mtags.find(_B_tag) == mtags.end())
+      {
         object->useMatrixTag(_A_tag);
+        // Noneigen Kernels
+        object->useMatrixTag(_precond_tag);
+      }
       else
         object->useVectorTag(_Bx_tag);
     }
@@ -173,10 +215,6 @@ NonlinearEigenSystem::turnOffOnEigenNodalBCs(bool on)
 void
 NonlinearEigenSystem::solve()
 {
-  // If we want to take into consideration eigen kernels, and then add the eigen tag.
-  if (_precond_matrix_includes_eigen)
-    _precond_tags.insert(_B_tag);
-
   // Clear the iteration counters
   _current_l_its.clear();
   _current_nl_its = 0;