Eigenvalue output and colorful residuals

Closes idaholab#14292

framework/include/outputs/OutputWarehouse.h

@@ -194,6 +194,12 @@ class OutputWarehouse
   /// Reset the output system
   void reset();
 
+  /**
+   * Calls the timestepSetup function for each of the output objects
+   * @see FEProblemBase::solve()
+   */
+  void solveSetup();
+
 private:
   /**
    * Calls the outputStep method for each output object
@@ -247,12 +253,6 @@ class OutputWarehouse
    */
   void timestepSetup();
 
-  /**
-   * Calls the timestepSetup function for each of the output objects
-   * @see FEProblemBase::solve()
-   */
-  void solveSetup();
-
   /**
    * Calls the jacobianSetup function for each of the output objects
    * @see FEProblemBase::computeJacobian

framework/include/problems/EigenProblem.h

@@ -138,6 +138,11 @@ class EigenProblem : public FEProblemBase
    * Which eigenvalue is active
    */
   unsigned int activeEigenvalueIndex() { return _active_eigen_index; }
+
+  const ConsoleStream & console() { return _console; }
+
+  virtual void initPetscOutput() override;
+
 #endif
 
 protected:

framework/include/problems/FEProblemBase.h

@@ -526,7 +526,7 @@ class FEProblemBase : public SubProblem, public Restartable
   /**
    * Reinitialize petsc output for proper linear/nonlinear iteration display
    */
-  void initPetscOutput();
+  virtual void initPetscOutput();
 
 #ifdef LIBMESH_HAVE_PETSC
   /**

framework/include/systems/DumpObjectsNonlinearSystem.h

@@ -31,6 +31,7 @@ class DumpObjectsNonlinearSystem : public NonlinearSystemBase
   virtual void stopSolve() override {}
   virtual bool converged() override { return true; }
   virtual NumericVector<Number> & RHS() override { return *_dummy; }
+  virtual SNES getSNES() override { return nullptr; }
 
   virtual unsigned int getCurrentNonlinearIterationNumber() override { return 0; }
   virtual void setupFiniteDifferencedPreconditioner() override {}

framework/include/systems/NonlinearEigenSystem.h

@@ -93,6 +93,8 @@ class NonlinearEigenSystem : public NonlinearSystemBase
 
   virtual NonlinearSolver<Number> * nonlinearSolver() override;
 
+  virtual SNES getSNES() override;
+
   virtual TransientEigenSystem & sys() { return _transient_sys; }
 
   /**

framework/include/systems/NonlinearSystem.h

@@ -64,6 +64,8 @@ class NonlinearSystem : public NonlinearSystemBase
     return _transient_sys.nonlinear_solver.get();
   }
 
+  virtual SNES getSNES() override;
+
   virtual TransientNonlinearImplicitSystem & sys() { return _transient_sys; }
 
   virtual void attachPreconditioner(Preconditioner<Number> * preconditioner) override;

framework/include/systems/NonlinearSystemBase.h

@@ -85,6 +85,8 @@ class NonlinearSystemBase : public SystemBase, public PerfGraphInterface
 
   virtual NonlinearSolver<Number> * nonlinearSolver() = 0;
 
+  virtual SNES getSNES() = 0;
+
   virtual unsigned int getCurrentNonlinearIterationNumber() = 0;
 
   /**

framework/include/utils/SlepcSupport.h

@@ -21,7 +21,7 @@
 /* We need this in order to implement moose PC */
 #include <petsc/private/pcimpl.h>
 #include <slepceps.h>
-#include <slepc/private/epsimpl.h>  
+#include <slepc/private/epsimpl.h>
 /* In order to use libMesh preconditioner */
 #include "libmesh/linear_solver.h"
 #include "libmesh/preconditioner.h"
@@ -45,7 +45,7 @@ void storeSlepcOptions(FEProblemBase & fe_problem, const InputParameters & param
 void storeSlepcEigenProblemOptions(EigenProblem & eigen_problem, const InputParameters & params);
 void slepcSetOptions(EigenProblem & eigen_problem, const InputParameters & params);
 void setSlepcEigenSolverTolerances(EigenProblem & eigen_problem, const InputParameters & params);
-void setSlepcOutputOptions(EigenProblem & eigen_problem);
+void setSlepcOutputOptions();
 void setFreeNonlinearPowerIterations(unsigned int free_power_iterations);
 void clearFreeNonlinearPowerIterations(const InputParameters & params);
 void moosePetscSNESFormMatrixTag(SNES snes, Vec x, Mat mat, void * ctx, TagID tag);
@@ -57,6 +57,15 @@ PetscErrorCode mooseSlepcEigenFormFunctionA(SNES snes, Vec x, Vec r, void * ctx)
 PetscErrorCode mooseSlepcEigenFormFunctionB(SNES snes, Vec x, Vec r, void * ctx);
 PetscErrorCode mooseSlepcEigenFormFunctionAB(SNES snes, Vec x, Vec Ax, Vec Bx, void * ctx);
 PetscErrorCode mooseSlepcStoppingTest(EPS eps,PetscInt its,PetscInt max_it,PetscInt nconv,PetscInt nev,EPSConvergedReason *reason,void *ctx);
+PetscErrorCode epsGetSNES(EPS eps, SNES * snes);
+PetscErrorCode mooseSlepcEPSMonitor(EPS eps,
+                                    int its,
+                                    int nconv,
+                                    PetscScalar * eigr,
+                                    PetscScalar * eigi,
+                                    PetscReal * errest,
+                                    int nest,
+                                    void * mctx);
 
 void attachCallbacksToMat(EigenProblem & eigen_problem, Mat mat, bool eigen);
 

framework/src/executioners/Eigenvalue.C

@@ -52,13 +52,12 @@ Eigenvalue::validParams()
       "normal_factor", 1.0, "Normalize eigenvector to make a defined norm equal to this factor");
 
   params.addParam<bool>("auto_initialization",
-                        false,
+                        true,
                         "If true, we will set an initial eigen vector in moose, otherwise EPS "
                         "solver will initial eigen vector");
 
-  params.addParam<bool>("newton_inverse_power",
-                        true,
-                        "If Newton and Inverse Power is combined in SLEPc side");
+  params.addParam<bool>(
+      "newton_inverse_power", false, "If Newton and Inverse Power is combined in SLEPc side");
 
 // Add slepc options and eigen problems
 #ifdef LIBMESH_HAVE_SLEPC
@@ -122,6 +121,8 @@ Eigenvalue::init()
       if (_eigen_problem.needInitializeEigenVector())
         _eigen_problem.initEigenvector(1.0);
 
+      _console << " Free power iteration starts" << std::endl;
+
       // Call solver
       _eigen_problem.solve();
       // Clear free power iterations
@@ -144,6 +145,9 @@ Eigenvalue::execute()
     _eigen_problem.doInitialFreePowerIteration(true);
     // Set free power iterations
     setFreeNonlinearPowerIterations(extra_power_iterations);
+
+    _console << " Extra Free power iteration starts" << std::endl;
+
     // Call solver
     _eigen_problem.solve();
     // Clear free power iterations
@@ -152,6 +156,11 @@ Eigenvalue::execute()
     _eigen_problem.doInitialFreePowerIteration(false);
   }
 
+  if (_eigen_problem.solverParams()._eigen_solve_type != Moose::EST_NONLINEAR_POWER)
+    _console << " Nonlinear Newton iteration starts" << std::endl;
+  else
+    _console << " Nonlinear power iteration starts" << std::endl;
+
   Steady::execute();
 }
 

framework/src/outputs/PetscOutput.C

@@ -134,9 +134,7 @@ PetscOutput::solveSetup()
 
   // Extract the non-linear and linear solvers from PETSc
   NonlinearSystemBase & nl = _problem_ptr->getNonlinearSystemBase();
-  PetscNonlinearSolver<Number> * petsc_solver =
-      dynamic_cast<PetscNonlinearSolver<Number> *>(nl.nonlinearSolver());
-  SNES snes = petsc_solver->snes();
+  SNES snes = nl.getSNES();
   KSP ksp;
   SNESGetKSP(snes, &ksp);
 

framework/src/problems/EigenProblem.C

@@ -411,3 +411,9 @@ EigenProblem::needInitializeEigenVector()
 {
   return _auto_initilize_eigen_vector && isNonlinearEigenvalueSolver();
 }
+
+void
+EigenProblem::initPetscOutput()
+{
+  _app.getOutputWarehouse().solveSetup();
+}

framework/src/systems/NonlinearEigenSystem.C

@@ -397,6 +397,26 @@ NonlinearEigenSystem::nonlinearSolver()
   return NULL;
 }
 
+SNES
+NonlinearEigenSystem::getSNES()
+{
+  SlepcEigenSolver<Number> * solver =
+      libmesh_cast_ptr<SlepcEigenSolver<Number> *>(&(*_transient_sys.eigen_solver));
+
+  if (!solver)
+    mooseError("There is no a eigen solver");
+
+  if (_eigen_problem.isNonlinearEigenvalueSolver())
+  {
+    EPS eps = solver->eps();
+    SNES snes = nullptr;
+    Moose::SlepcSupport::epsGetSNES(eps, &snes);
+    return snes;
+  }
+  else
+    mooseError("There is no SNES in linear eigen solver");
+}
+
 void
 NonlinearEigenSystem::checkIntegrity()
 {

framework/src/systems/NonlinearSystem.C

@@ -446,3 +446,15 @@ NonlinearSystem::computeScalingResidual()
 {
   _fe_problem.computeResidualSys(_transient_sys, *_current_solution, RHS());
 }
+
+SNES
+NonlinearSystem::getSNES()
+{
+  PetscNonlinearSolver<Number> * petsc_solver =
+      dynamic_cast<PetscNonlinearSolver<Number> *>(nonlinearSolver());
+
+  if (petsc_solver)
+    return petsc_solver->snes();
+  else
+    mooseError("It is not a petsc nonlinear solver");
+}

framework/src/utils/SlepcEigenSolverConfiguration.C

@@ -26,10 +26,20 @@ SlepcEigenSolverConfiguration::SlepcEigenSolverConfiguration(EigenProblem & eige
 void
 SlepcEigenSolverConfiguration::configure_solver()
 {
-  PetscErrorCode ierr;
-
-  std::cout<<" configure_solver "<<std::endl;
-
-  ierr = EPSSetStoppingTestFunction(_slepc_solver.eps(),Moose::SlepcSupport::mooseSlepcStoppingTest,&_eigen_problem,NULL);
-  LIBMESH_CHKERR(ierr);
+  PetscErrorCode ierr = 0;
+
+  if (_eigen_problem.isNonlinearEigenvalueSolver())
+  {
+    _eigen_problem.initPetscOutput();
+
+    ierr = EPSSetStoppingTestFunction(
+        _slepc_solver.eps(), Moose::SlepcSupport::mooseSlepcStoppingTest, &_eigen_problem, NULL);
+    LIBMESH_CHKERR(ierr);
+
+    ierr = EPSMonitorCancel(_slepc_solver.eps());
+    LIBMESH_CHKERR(ierr);
+    ierr = EPSMonitorSet(
+        _slepc_solver.eps(), Moose::SlepcSupport::mooseSlepcEPSMonitor, &_eigen_problem, NULL);
+    LIBMESH_CHKERR(ierr);
+  }
 }

framework/src/utils/SlepcSupport.C

@@ -269,52 +269,10 @@ setEigenProblemOptions(SolverParams & solver_params)
 }
 
 void
-setNewtonOutputOptions()
-{
-  Moose::PetscSupport::setSinglePetscOption("-init_eps_monitor_conv");
-  Moose::PetscSupport::setSinglePetscOption("-init_eps_monitor");
-  Moose::PetscSupport::setSinglePetscOption("-eps_power_snes_monitor");
-  Moose::PetscSupport::setSinglePetscOption("-eps_power_ksp_monitor");
-  Moose::PetscSupport::setSinglePetscOption("-init_eps_power_snes_monitor");
-  Moose::PetscSupport::setSinglePetscOption("-init_eps_power_ksp_monitor");
-}
-
-void
-setSlepcOutputOptions(EigenProblem & eigen_problem)
+setSlepcOutputOptions()
 {
   Moose::PetscSupport::setSinglePetscOption("-eps_monitor_conv");
   Moose::PetscSupport::setSinglePetscOption("-eps_monitor");
-  switch (eigen_problem.solverParams()._eigen_solve_type)
-  {
-    case Moose::EST_NONLINEAR_POWER:
-      Moose::PetscSupport::setSinglePetscOption("-eps_power_snes_monitor");
-      Moose::PetscSupport::setSinglePetscOption("-eps_power_ksp_monitor");
-      break;
-
-    case Moose::EST_NEWTON:
-      setNewtonOutputOptions();
-      break;
-
-    case Moose::EST_PJFNK:
-      setNewtonOutputOptions();
-      break;
-
-    case Moose::EST_JFNK:
-      setNewtonOutputOptions();
-      break;
-
-    case Moose::EST_POWER:
-      break;
-
-    case Moose::EST_ARNOLDI:
-      break;
-
-    case Moose::EST_KRYLOVSCHUR:
-      break;
-
-    case Moose::EST_JACOBI_DAVIDSON:
-      break;
-  }
 }
 
 void
@@ -504,7 +462,7 @@ slepcSetOptions(EigenProblem & eigen_problem, const InputParameters & params)
   setEigenProblemOptions(eigen_problem.solverParams());
   setWhichEigenPairsOptions(eigen_problem.solverParams());
   setSlepcEigenSolverTolerances(eigen_problem, params);
-  setSlepcOutputOptions(eigen_problem);
+  setSlepcOutputOptions();
   Moose::PetscSupport::addPetscOptionsFromCommandline();
 }
 
@@ -1014,9 +972,48 @@ mooseSlepcStoppingTest(EPS eps,PetscInt its,PetscInt max_it,PetscInt nconv,Petsc
     *reason = EPS_CONVERGED_USER;
     eps->nconv = 1;
   }
+  return 0;
+}
+
+PetscErrorCode
+epsGetSNES(EPS eps, SNES * snes)
+{
+  PetscErrorCode ierr;
+  PetscBool same, nonlinear;
+
+  ierr = PetscObjectTypeCompare((PetscObject)eps, EPSPOWER, &same);
+  LIBMESH_CHKERR(ierr);
+
+  if (!same)
+    mooseError("It is not eps power, and there is no snes");
+
+  ierr = EPSPowerGetNonlinear(eps, &nonlinear);
+  LIBMESH_CHKERR(ierr);
+
+  if (!nonlinear)
+    mooseError("It is not a nonlinear eigen solver");
+
+  ierr = EPSPowerGetSNES(eps, snes);
+  LIBMESH_CHKERR(ierr);
+
+  return 0;
+}
+
+PetscErrorCode
+mooseSlepcEPSMonitor(EPS /*eps*/,
+                     int its,
+                     int /*nconv*/,
+                     PetscScalar * eigr,
+                     PetscScalar * /*eigi*/,
+                     PetscReal * /*errest*/,
+                     int /*nest*/,
+                     void * mctx)
+{
 
+  EigenProblem * eigen_problem = static_cast<EigenProblem *>(mctx);
+  auto & console = eigen_problem->console();
 
-  std::cout<<"its "<<its<<" max_it "<<max_it<<" nconv "<<nconv<<" nev "<<nev<<" reason "<<*reason<<std::endl;
+  console << " Iteration " << its << std::setprecision(8) << " eigenvalue = " << *eigr << std::endl;
 
   return 0;
 }

test/tests/problems/eigen_problem/arraykernels/tests

@@ -26,7 +26,7 @@
     csvdiff = 'ne_two_variables_out_eigenvalues_0001.csv'
     cli_args = 'Executioner/precond_matrix_includes_eigen=true -eps_power_pc_type  lu -init_eps_power_pc_type lu -eps_power_pc_factor_mat_solver_type superlu_dist'
     prereq = 'two_variables'
-    expect_out = '6 KSP Residual norm'
+    expect_out = '6 Linear |R|'
     slepc = true
     slepc_version = '>=3.11.0'
     requirement = "Eigenvalue system should be able to handle multiple variables"

test/tests/problems/eigen_problem/initial_condition/ne_ic.i

@@ -51,10 +51,7 @@
 
 [Executioner]
   type = Eigenvalue
-  matrix_free = true
-  solve_type = NEWTON
-  auto_initialization = true
-  newton_inverse_power = false
+  solve_type = PJFNK
   eigen_problem_type = GEN_NON_HERMITIAN
 []
 

test/tests/problems/eigen_problem/preconditioners/tests

@@ -49,7 +49,7 @@
     csvdiff = 'ne_pbp_out_eigenvalues_0001.csv'
     cli_args = 'Executioner/precond_matrix_free=true Executioner/precond_matrix_includes_eigen=true'
     prereq = 'newton_pbp_shell_precond'
-    expect_out = '6 KSP Residual norm'
+    expect_out = '6 Linear |R|'
     slepc_version = '>=3.8.0'
     requirement = "Eigen solver should work with a physics-based preconditioner with including eigen kernels in the preconditioning matrix"
   [../]