Merge pull request #25204 from lindsayad/slepc-normalization

custom eigenvalue normalization for slepc

framework/include/kernels/MassEigenKernel.h

@@ -21,4 +21,5 @@ class MassEigenKernel : public EigenKernel
 protected:
   virtual Real computeQpResidual() override;
   virtual Real computeQpJacobian() override;
+  const Real _coef;
 };

framework/include/problems/EigenProblem.h

@@ -224,21 +224,15 @@ class EigenProblem : public FEProblemBase
    */
   void wereMatricesFormed(bool mf) { _matrices_formed = mf; }
 
-private:
   /**
-   * Do some free/extra power iterations
+   * Form the Bx norm
    */
-  void doFreeNonlinearPowerIterations(unsigned int free_power_iterations);
+  Real formNorm();
 
   /**
-   * Adjust eigen vector by either scaling the existing values or setting new values
-   * The operations are applied for only eigen variables
+   * Whether a Bx norm postprocessor has been provided
    */
-  void adjustEigenVector(const Real value, bool scaling);
-
-#endif
-
-  using FEProblemBase::_nl;
+  bool bxNormProvided() const { return _bx_norm_name.has_value(); }
 
 protected:
   unsigned int _n_eigen_pairs_required;
@@ -275,6 +269,26 @@ class EigenProblem : public FEProblemBase
   bool & _first_solve;
   /// A value used for initial normalization
   Real _initial_eigenvalue;
+
+private:
+  /**
+   * Do some free/extra power iterations
+   */
+  void doFreeNonlinearPowerIterations(unsigned int free_power_iterations);
+
+  /**
+   * Adjust eigen vector by either scaling the existing values or setting new values
+   * The operations are applied for only eigen variables
+   */
+  void adjustEigenVector(const Real value, bool scaling);
+
+  /// The name of the Postprocessor providing the Bx norm. This may be empty in which case the
+  /// default L2 norm of Bx will be used as the Bx norm
+  const std::optional<PostprocessorName> _bx_norm_name;
+
+#endif
+
+  using FEProblemBase::_nl;
 };
 
 #ifdef LIBMESH_HAVE_SLEPC

framework/src/kernels/MassEigenKernel.C

@@ -15,21 +15,26 @@ InputParameters
 MassEigenKernel::validParams()
 {
   InputParameters params = EigenKernel::validParams();
-  params.addClassDescription("An eigenkernel with weak form $\\lambda(\\psi_i, -u_h)$ where "
-                             "$\\lambda$ is the eigenvalue.");
+  params.addClassDescription(
+      "An eigenkernel with weak form $\\lambda(\\psi_i, -u_h \\times coeff)$ where "
+      "$\\lambda$ is the eigenvalue.");
+  params.addParam<Real>("coefficient", 1.0, "Coefficient multiplying the term");
   return params;
 }
 
-MassEigenKernel::MassEigenKernel(const InputParameters & parameters) : EigenKernel(parameters) {}
+MassEigenKernel::MassEigenKernel(const InputParameters & parameters)
+  : EigenKernel(parameters), _coef(this->template getParam<Real>("coefficient"))
+{
+}
 
 Real
 MassEigenKernel::computeQpResidual()
 {
-  return -_u[_qp] * _test[_i][_qp];
+  return -_coef * _u[_qp] * _test[_i][_qp];
 }
 
 Real
 MassEigenKernel::computeQpJacobian()
 {
-  return -_phi[_j][_qp] * _test[_i][_qp];
+  return -_coef * _phi[_j][_qp] * _test[_i][_qp];
 }

framework/src/problems/EigenProblem.C

@@ -42,6 +42,7 @@ EigenProblem::validParams()
       "active_eigen_index",
       0,
       "Which eigenvector is used to compute residual and also associated to nonlinear variable");
+  params.addParam<PostprocessorName>("bx_norm", "A postprocessor describing the norm of Bx");
 
   return params;
 }
@@ -60,7 +61,10 @@ EigenProblem::EigenProblem(const InputParameters & parameters)
     _constant_matrices(false),
     _has_normalization(false),
     _normal_factor(1.0),
-    _first_solve(declareRestartableData<bool>("first_solve", true))
+    _first_solve(declareRestartableData<bool>("first_solve", true)),
+    _bx_norm_name(isParamValid("bx_norm")
+                      ? std::make_optional(getParam<PostprocessorName>("bx_norm"))
+                      : std::nullopt)
 {
 #ifdef LIBMESH_HAVE_SLEPC
   if (_nl_sys_names.size() > 1)
@@ -416,7 +420,7 @@ EigenProblem::preScaleEigenVector(const std::pair<Real, Real> & eig)
   // Eigenvalue magnitude
   Real v = std::sqrt(eig.first * eig.first + eig.second * eig.second);
   // Scaling factor
-  Real factor = 1 / v / _nl_eigen->residualVectorBX().l2_norm();
+  Real factor = 1 / v / (bxNormProvided() ? formNorm() : _nl_eigen->residualVectorBX().l2_norm());
   // Scale eigenvector
   if (!MooseUtils::absoluteFuzzyEqual(factor, 1))
     scaleEigenvector(factor);
@@ -637,4 +641,11 @@ EigenProblem::initPetscOutput()
   _app.getOutputWarehouse().solveSetup();
 }
 
+Real
+EigenProblem::formNorm()
+{
+  mooseAssert(_bx_norm_name,
+              "We should not get here unless a bx_norm postprocessor has been provided");
+  return getPostprocessorValueByName(*_bx_norm_name);
+}
 #endif

framework/src/utils/SlepcSupport.C

@@ -356,6 +356,9 @@ setFreeNonlinearPowerIterations(unsigned int free_power_iterations)
   // -snes_no_convergence_test is a perfect option, but it was removed from PETSc
   Moose::PetscSupport::setSinglePetscOption("-snes_rtol", "0.99999999999");
   Moose::PetscSupport::setSinglePetscOption("-eps_max_it", stringify(free_power_iterations));
+  // We always want the number of free power iterations respected so we don't want to stop early if
+  // we've satisfied a convergence criterion. Consequently we make this tolerance very tight
+  Moose::PetscSupport::setSinglePetscOption("-eps_tol", "1e-50");
 }
 
 void
@@ -367,6 +370,7 @@ clearFreeNonlinearPowerIterations(const InputParameters & params)
                                             stringify(params.get<unsigned int>("nl_max_its")));
   Moose::PetscSupport::setSinglePetscOption("-snes_rtol",
                                             stringify(params.get<Real>("nl_rel_tol")));
+  Moose::PetscSupport::setSinglePetscOption("-eps_tol", stringify(params.get<Real>("eigen_tol")));
 }
 
 void
@@ -962,28 +966,60 @@ mooseSlepcEigenFormFunctionAB(SNES /*snes*/, Vec x, Vec Ax, Vec Bx, void * ctx)
   PetscFunctionReturn(0);
 }
 
+PetscErrorCode
+mooseSlepcEigenFormNorm(SNES /*snes*/, Vec /*Bx*/, PetscReal * norm, void * ctx)
+{
+  PetscFunctionBegin;
+  auto * const eigen_problem = static_cast<EigenProblem *>(ctx);
+  *norm = eigen_problem->formNorm();
+  PetscFunctionReturn(0);
+}
+
 void
 attachCallbacksToMat(EigenProblem & eigen_problem, Mat mat, bool eigen)
 {
+  // Recall that we are solving the potentially nonlinear problem:
+  // F(x) = A(x) - \lambda B(x) = 0
+  //
+  // To solve this, we can use Newton's method: J \Delta x = -F
+  // Generally we will approximate J using matrix free methods. However, in order to solve the
+  // linearized system efficiently, we typically will need preconditioning. Typically we will build
+  // the preconditioner only from A, but we also have the option to include information from B
+
+  // Attach the Jacobian computation function. If \p mat is the "eigen" matrix corresponding to B,
+  // then attach our JacobianB computation routine, else the matrix corresponds to A, and we attach
+  // the JacobianA computation routine
   PetscObjectComposeFunction((PetscObject)mat,
                              "formJacobian",
                              eigen ? Moose::SlepcSupport::mooseSlepcEigenFormJacobianB
                                    : Moose::SlepcSupport::mooseSlepcEigenFormJacobianA);
+
+  // Attach the residual computation function. If \p mat is the "eigen" matrix corresponding to B,
+  // then attach our FunctionB computation routine, else the matrix corresponds to A, and we attach
+  // the FunctionA computation routine
   PetscObjectComposeFunction((PetscObject)mat,
                              "formFunction",
                              eigen ? Moose::SlepcSupport::mooseSlepcEigenFormFunctionB
                                    : Moose::SlepcSupport::mooseSlepcEigenFormFunctionA);
 
+  // It's also beneficial to be able to evaluate both A and B residuals at once
   PetscObjectComposeFunction(
       (PetscObject)mat, "formFunctionAB", Moose::SlepcSupport::mooseSlepcEigenFormFunctionAB);
 
+  // Users may choose to provide a custom measure of the norm of B (Bx for a linear system)
+  if (eigen_problem.bxNormProvided())
+    PetscObjectComposeFunction(
+        (PetscObject)mat, "formNorm", Moose::SlepcSupport::mooseSlepcEigenFormNorm);
+
+  // Finally we need to attach the "context" object, which is our EigenProblem, to the matrices so
+  // that eventually when we get callbacks from SLEPc we can call methods on the EigenProblem
   PetscContainer container;
   PetscContainerCreate(eigen_problem.comm().get(), &container);
   PetscContainerSetPointer(container, &eigen_problem);
-  PetscObjectCompose((PetscObject)mat, "formJacobianCtx", nullptr);
   PetscObjectCompose((PetscObject)mat, "formJacobianCtx", (PetscObject)container);
-  PetscObjectCompose((PetscObject)mat, "formFunctionCtx", nullptr);
   PetscObjectCompose((PetscObject)mat, "formFunctionCtx", (PetscObject)container);
+  if (eigen_problem.bxNormProvided())
+    PetscObjectCompose((PetscObject)mat, "formNormCtx", (PetscObject)container);
   PetscContainerDestroy(&container);
 }
 
@@ -1086,7 +1122,8 @@ PCCreate_MoosePC(PC pc)
   PetscFunctionReturn(0);
 }
 
-PetscErrorCode PCDestroy_MoosePC(PC /*pc*/)
+PetscErrorCode
+PCDestroy_MoosePC(PC /*pc*/)
 {
   PetscFunctionBegin;
   /* We do not need to do anything right now, but later we may have some data we need to free here

test/tests/executioners/eigen_convergence/a.i

@@ -0,0 +1,84 @@
+[Mesh]
+  [gen]
+    type = GeneratedMeshGenerator
+    dim = 2
+    xmin = 0
+    xmax = 160
+    ymin = 0
+    ymax = 160
+    nx = 8
+    ny = 8
+  []
+  uniform_refine = 0
+[]
+
+[Variables]
+  [u]
+  []
+[]
+
+[Kernels]
+  [diffusion]
+    type = MatDiffusion
+    diffusivity = diffusivity
+    variable = u
+  []
+  [reaction]
+    type = CoefReaction
+    coefficient = 0.01
+    variable = u
+  []
+  [rhs]
+    type = CoefReaction
+    extra_vector_tags = 'eigen'
+    coefficient = -0.01
+    variable = u
+  []
+[]
+
+[BCs]
+  [robin]
+    type = VacuumBC
+    boundary = 'left bottom'
+    variable = u
+  []
+[]
+
+[Materials]
+  [nm]
+    type = GenericConstantMaterial
+    block = 0
+    prop_names = 'diffusivity'
+    prop_values = 0.333333333333333333
+  []
+[]
+
+[VectorPostprocessors]
+  [eigen]
+    type = Eigenvalues
+    inverse_eigenvalue = true
+  []
+[]
+
+[Postprocessors]
+  [fluxintegral]
+    type = ElementIntegralVariablePostprocessor
+    variable = u
+    execute_on = linear
+  []
+[]
+
+[Problem]
+  type = EigenProblem
+[]
+
+[Executioner]
+  type = Eigenvalue
+  solve_type = PJFNK
+  free_power_iterations = 4
+  nl_abs_tol = 2e-10
+[]
+
+[Outputs]
+  csv = true
+[]

test/tests/executioners/eigen_convergence/b.i

@@ -0,0 +1,73 @@
+[Mesh]
+  [gen]
+    type = GeneratedMeshGenerator
+    dim = 2
+    xmin = 0
+    xmax = 160
+    ymin = 0
+    ymax = 160
+    nx = 8
+    ny = 8
+  []
+  uniform_refine = 0
+[]
+
+[Variables]
+  [u]
+  []
+[]
+
+[Kernels]
+  [diffusion]
+    type = MatDiffusion
+    diffusivity = diffusivity
+    variable = u
+  []
+  [reaction]
+    type = CoefReaction
+    coefficient = 0.01
+    variable = u
+  []
+  [rhs]
+    type = MassEigenKernel
+    coefficient = 0.01
+    variable = u
+  []
+[]
+
+[BCs]
+  [robin]
+    type = VacuumBC
+    boundary = 'left bottom'
+    variable = u
+  []
+[]
+
+[Materials]
+  [nm]
+    type = GenericConstantMaterial
+    block = 0
+    prop_names = 'diffusivity'
+    prop_values = 0.333333333333333333
+  []
+[]
+
+[Postprocessors]
+  [fluxintegral]
+    type = ElementIntegralVariablePostprocessor
+    variable = u
+    execute_on = linear
+  []
+[]
+
+[Executioner]
+  type = NonlinearEigen
+  bx_norm = fluxintegral
+  solve_type = PJFNK
+  free_power_iterations = 4
+  nl_abs_tol = 2e-10
+[]
+
+[Outputs]
+  csv = true
+[]

test/tests/executioners/eigen_convergence/gold/a_out_eigen_0001.csv

@@ -0,0 +1,2 @@
+eigen_values_imag,eigen_values_real
+0,0.99364796564363

test/tests/executioners/eigen_convergence/gold/b_out.csv

@@ -0,0 +1,4 @@
+time,fluxintegral,Executioner/eigenvalue
+0,1,1
+1,0.97869657767363,0.97869657767363
+2,0.99364796612084,0.99364796612084