apply bx norm factor for eigen kernel tagging idaholab#28984

framework/include/interfaces/TaggingInterface.h

@@ -27,6 +27,7 @@ class MooseObject;
 class SubProblem;
 class Assembly;
 class ReferenceResidualProblem;
+class EigenProblem;
 
 template <typename T>
 InputParameters validParams();
@@ -395,6 +396,11 @@ class TaggingInterface
   /// this data member to avoid constant dynamic heap allocations
   std::vector<Real> _absolute_residuals;
 
+  /// Pointer to the eigen problem for an eigen object
+  const EigenProblem * _eigen_problem = nullptr;
+  /// Pointer to the postprocessor for the Bx norm
+  const PostprocessorValue * _bx_norm = nullptr;
+
   friend void NonlinearSystemBase::constraintJacobians(bool);
 };
 

framework/include/problems/EigenProblem.h

@@ -49,6 +49,7 @@ class EigenProblem : public FEProblemBase
 
   NonlinearEigenSystem & getNonlinearEigenSystem(const unsigned int nl_sys_num);
   NonlinearEigenSystem & getCurrentNonlinearEigenSystem();
+  const NonlinearEigenSystem & getCurrentNonlinearEigenSystem() const;
 
   virtual void checkProblemIntegrity() override;
 
@@ -227,7 +228,7 @@ class EigenProblem : public FEProblemBase
   /**
    * Form the Bx norm
    */
-  Real formNorm();
+  const Real & formNorm() const;
 
   /**
    * Whether a Bx norm postprocessor has been provided
@@ -239,6 +240,8 @@ class EigenProblem : public FEProblemBase
    */
   void setBxNorm(const PostprocessorName & bx_norm) { _bx_norm_name = bx_norm; }
 
+  const PostprocessorName & bxNormName() const { return _bx_norm_name.value(); }
+
 protected:
   unsigned int _n_eigen_pairs_required;
   bool _generalized_eigenvalue_problem;
@@ -291,6 +294,8 @@ class EigenProblem : public FEProblemBase
   /// default L2 norm of Bx will be used as the Bx norm
   std::optional<PostprocessorName> _bx_norm_name;
 
+  /// Pointer to the postprocessor providing the Bx norm
+  const PostprocessorValue * _bx_norm = nullptr;
 #endif
 
   using FEProblemBase::_nl;
@@ -312,4 +317,10 @@ EigenProblem::getCurrentNonlinearEigenSystem()
   return *_nl_eigen;
 }
 
+inline const NonlinearEigenSystem &
+EigenProblem::getCurrentNonlinearEigenSystem() const
+{
+  return *_nl_eigen;
+}
+
 #endif

framework/src/interfaces/TaggingInterface.C

@@ -12,6 +12,8 @@
 #include "FEProblem.h"
 #include "Assembly.h"
 #include "ReferenceResidualProblem.h"
+#include "EigenProblem.h"
+#include "NonlinearEigenSystem.h"
 
 #include "libmesh/dense_vector.h"
 
@@ -147,6 +149,19 @@ TaggingInterface::TaggingInterface(const MooseObject * moose_object)
     _non_ref_vector_tags = _vector_tags;
     _non_ref_abs_vector_tags = _abs_vector_tags;
   }
+
+  // we need to scale the eigen kernel residual when tagging to a vector different than Bx vector
+  if (const auto * const eigen_problem = dynamic_cast<const EigenProblem *>(fe_problem))
+  {
+    const auto & eigen_system = eigen_problem->getCurrentNonlinearEigenSystem();
+    if (_vector_tags.count(eigen_system.eigenVectorTag()))
+    {
+      _eigen_problem = eigen_problem;
+
+      if (_eigen_problem->bxNormProvided())
+        _bx_norm = &_eigen_problem->getPostprocessorValueByName(eigen_problem->bxNormName());
+    }
+  }
 }
 
 void
@@ -359,8 +374,15 @@ TaggingInterface::prepareMatrixTagLower(Assembly & assembly,
 void
 TaggingInterface::accumulateTaggedLocalResidual()
 {
+  if (_eigen_problem)
+  {
+    Real fac = (_eigen_problem->negativeSignEigenKernel() ? 1 : -1) * (_bx_norm ? *_bx_norm : 1);
+    _local_re.scale(1 / fac);
+  }
+
   for (auto & re : _re_blocks)
     *re += _local_re;
+
   for (auto & absre : _absre_blocks)
     for (const auto i : index_range(_local_re))
       (*absre)(i) += std::abs(_local_re(i));
@@ -369,8 +391,15 @@ TaggingInterface::accumulateTaggedLocalResidual()
 void
 TaggingInterface::assignTaggedLocalResidual()
 {
+  if (_eigen_problem)
+  {
+    Real fac = (_eigen_problem->negativeSignEigenKernel() ? 1 : -1) * (_bx_norm ? *_bx_norm : 1);
+    _local_re.scale(1 / fac);
+  }
+
   for (auto & re : _re_blocks)
     *re = _local_re;
+
   for (auto & absre : _absre_blocks)
     for (const auto i : index_range(_local_re))
       (*absre)(i) = std::abs(_local_re(i));
@@ -379,8 +408,12 @@ TaggingInterface::assignTaggedLocalResidual()
 void
 TaggingInterface::accumulateTaggedLocalMatrix()
 {
-  for (auto & ke : _ke_blocks)
-    *ke += _local_ke;
+  if (_eigen_problem && _eigen_problem->negativeSignEigenKernel())
+    for (auto & ke : _ke_blocks)
+      *ke -= _local_ke;
+  else
+    for (auto & ke : _ke_blocks)
+      *ke += _local_ke;
 }
 
 void
@@ -396,8 +429,12 @@ TaggingInterface::accumulateTaggedLocalMatrix(Assembly & assembly,
     _ke_blocks[i] = &assembly.jacobianBlock(ivar, jvar, Assembly::LocalDataKey{}, *mat_vector);
   mooseAssert(_ke_blocks[0]->m() == k.m() && _ke_blocks[0]->n() == k.n(),
               "Passed-in k must match the blocks we are about to sum into");
-  for (auto & ke : _ke_blocks)
-    *ke += k;
+  if (_eigen_problem && _eigen_problem->negativeSignEigenKernel())
+    for (auto & ke : _ke_blocks)
+      *ke -= k;
+  else
+    for (auto & ke : _ke_blocks)
+      *ke += k;
 }
 
 void
@@ -415,22 +452,37 @@ TaggingInterface::accumulateTaggedLocalMatrix(Assembly & assembly,
         &assembly.jacobianBlockNeighbor(type, ivar, jvar, Assembly::LocalDataKey{}, *mat_vector);
   mooseAssert(_ke_blocks[0]->m() == k.m() && _ke_blocks[0]->n() == k.n(),
               "Passed-in k must match the blocks we are about to sum into");
-  for (auto & ke : _ke_blocks)
-    *ke += k;
+  if (_eigen_problem && _eigen_problem->negativeSignEigenKernel())
+    for (auto & ke : _ke_blocks)
+      *ke -= k;
+  else
+    for (auto & ke : _ke_blocks)
+      *ke += k;
 }
 
 void
 TaggingInterface::accumulateTaggedNonlocalMatrix()
 {
-  for (auto & ke : _ke_blocks)
-    *ke += _nonlocal_ke;
+  if (_eigen_problem && _eigen_problem->negativeSignEigenKernel())
+    for (auto & ke : _ke_blocks)
+      *ke -= _nonlocal_ke;
+  else
+    for (auto & ke : _ke_blocks)
+      *ke += _nonlocal_ke;
 }
 
 void
 TaggingInterface::assignTaggedLocalMatrix()
 {
-  for (auto & ke : _ke_blocks)
-    *ke = _local_ke;
+  if (_eigen_problem && _eigen_problem->negativeSignEigenKernel())
+    for (auto & ke : _ke_blocks)
+    {
+      *ke = _local_ke;
+      ke->scale(-1);
+    }
+  else
+    for (auto & ke : _ke_blocks)
+      *ke = _local_ke;
 }
 
 TaggingInterface::~TaggingInterface() {}

framework/src/problems/EigenProblem.C

@@ -312,6 +312,14 @@ EigenProblem::computeResidualTag(const NumericVector<Number> & soln,
   setCurrentNonlinearSystem(_nl_eigen->number());
   computeResidualTags(_fe_vector_tags);
 
+  if (tag == _nl_eigen->eigenVectorTag())
+  {
+    if (bxNormProvided())
+      residual *= (-*_bx_norm);
+    else
+      residual *= -1;
+  }
+
   _nl_eigen->disassociateVectorFromTag(residual, tag);
 }
 
@@ -346,6 +354,21 @@ EigenProblem::computeResidualAB(const NumericVector<Number> & soln,
 
   computeResidualTags(_fe_vector_tags);
 
+  if (bxNormProvided())
+  {
+    if (tagA == _nl_eigen->eigenVectorTag())
+      residualA *= (-*_bx_norm);
+    if (tagB == _nl_eigen->eigenVectorTag())
+      residualB *= (-*_bx_norm);
+  }
+  else
+  {
+    if (tagA == _nl_eigen->eigenVectorTag())
+      residualA *= -1;
+    if (tagB == _nl_eigen->eigenVectorTag())
+      residualB *= -1;
+  }
+
   _nl_eigen->disassociateVectorFromTag(residualA, tagA);
   _nl_eigen->disassociateVectorFromTag(residualB, tagB);
 }
@@ -368,10 +391,7 @@ EigenProblem::computeResidualL2Norm()
   _nl_eigen->residualVectorBX() *= eigenvalue;
 
   // Compute entire residual
-  if (_negative_sign_eigen_kernel)
-    _nl_eigen->residualVectorAX() += _nl_eigen->residualVectorBX();
-  else
-    _nl_eigen->residualVectorAX() -= _nl_eigen->residualVectorBX();
+  _nl_eigen->residualVectorAX() -= _nl_eigen->residualVectorBX();
 
   return _nl_eigen->residualVectorAX().l2_norm();
 }
@@ -642,6 +662,9 @@ EigenProblem::init()
   _nl_eigen->sys().use_shell_precond_matrix(solverParams()._precond_matrix_free);
 #endif
 
+  if (bxNormProvided())
+    _bx_norm = &getPostprocessorValueByName(*_bx_norm_name);
+
   FEProblemBase::init();
 }
 
@@ -670,11 +693,11 @@ EigenProblem::initPetscOutputAndSomeSolverSettings()
   _app.getOutputWarehouse().solveSetup();
 }
 
-Real
-EigenProblem::formNorm()
+const Real &
+EigenProblem::formNorm() const
 {
   mooseAssert(_bx_norm_name,
               "We should not get here unless a bx_norm postprocessor has been provided");
-  return getPostprocessorValueByName(*_bx_norm_name);
+  return *_bx_norm;
 }
 #endif

framework/src/systems/NonlinearEigenSystem.C

@@ -81,14 +81,6 @@ assemble_matrix(EquationSystems & es, const std::string & system_name)
                          eigen_system.get_matrix_B(),
                          eigen_nl.nonEigenMatrixTag(),
                          eigen_nl.eigenMatrixTag());
-#if LIBMESH_HAVE_SLEPC
-    if (p->negativeSignEigenKernel())
-    {
-      auto ierr =
-          MatScale(static_cast<PetscMatrix<Number> &>(eigen_system.get_matrix_B()).mat(), -1.0);
-      LIBMESH_CHKERR(ierr);
-    }
-#endif
     return;
   }
 

framework/src/utils/SlepcSupport.C

@@ -866,12 +866,6 @@ mooseSlepcEigenFormJacobianB(SNES snes, Vec x, Mat jac, Mat pc, void * ctx)
 
   moosePetscSNESFormMatrixTag(snes, x, pc, sys.get_matrix_B(), ctx, eigen_nl.eigenMatrixTag());
 
-  if (eigen_problem->negativeSignEigenKernel())
-  {
-    ierr = MatScale(pc, -1.);
-    CHKERRQ(ierr);
-  }
-
   PetscFunctionReturn(PETSC_SUCCESS);
 }
 
@@ -955,12 +949,6 @@ mooseSlepcEigenFormFunctionB(SNES snes, Vec x, Vec r, void * ctx)
   else
     moosePetscSNESFormFunction(snes, x, r, ctx, eigen_nl.eigenVectorTag());
 
-  if (eigen_problem->negativeSignEigenKernel())
-  {
-    ierr = VecScale(r, -1.);
-    CHKERRQ(ierr);
-  }
-
   PetscFunctionReturn(PETSC_SUCCESS);
 }
 
@@ -1000,12 +988,6 @@ mooseSlepcEigenFormFunctionAB(SNES /*snes*/, Vec x, Vec Ax, Vec Bx, void * ctx)
     ierr = MatMult(B, x, Bx);
     CHKERRQ(ierr);
 
-    if (eigen_problem->negativeSignEigenKernel())
-    {
-      ierr = VecScale(Bx, -1.);
-      CHKERRQ(ierr);
-    }
-
     PetscFunctionReturn(PETSC_SUCCESS);
   }
 
@@ -1030,12 +1012,6 @@ mooseSlepcEigenFormFunctionAB(SNES /*snes*/, Vec x, Vec Ax, Vec Bx, void * ctx)
     sys.copy_super_to_sub(*BX, sub_Bx);
   }
 
-  if (eigen_problem->negativeSignEigenKernel())
-  {
-    ierr = VecScale(Bx, -1.);
-    CHKERRQ(ierr);
-  }
-
   PetscFunctionReturn(PETSC_SUCCESS);
 }
 
@@ -1125,12 +1101,6 @@ mooseMatMult_Eigen(Mat mat, Vec x, Vec r)
 
   evaluateResidual(*eigen_problem, x, r, eigen_nl.eigenVectorTag());
 
-  if (eigen_problem->negativeSignEigenKernel())
-  {
-    ierr = VecScale(r, -1.);
-    CHKERRQ(ierr);
-  }
-
   PetscFunctionReturn(PETSC_SUCCESS);
 }
 

test/tests/problems/eigen_problem/arraykernels/tests

@@ -29,7 +29,7 @@
       csvdiff = 'ne_two_variables_out_eigenvalues_0001.csv'
       cli_args = 'Executioner/precond_matrix_includes_eigen=true -pc_type  lu  lu '
                  '-pc_factor_mat_solver_type superlu_dist'
-      expect_out = '6 Linear \|R\|'
+      expect_out = '4 Linear \|R\|'
       slepc = true
       # Include B into the preconditioning matrix in general is a bad idea
       # If there is no good guess, free power iteration may converge a different mode