An attempt to add picard to eigenvalue system

Closes #12767
idaholab · Oct 17, 2019 · 9efa420 · 9efa420
1 parent 97326a7
commit 9efa420
Show file tree

Hide file tree

Showing 10 changed files with 197 additions and 8 deletions.
diff --git a/framework/include/problems/EigenProblem.h b/framework/include/problems/EigenProblem.h
@@ -11,6 +11,7 @@
 
 // MOOSE Includes
 #include "FEProblemBase.h"
+#include "Eigenvalue.h"
 
 // Forward declarations
 class EigenProblem;
@@ -55,6 +56,8 @@ class EigenProblem : public FEProblemBase
 #if LIBMESH_HAVE_SLEPC
   void setEigenproblemType(Moose::EigenProblemType eigen_problem_type);
 
+  virtual Real computeResidualL2Norm() override;
+
   /**
    * Form a Jacobian matrix for all kernels and BCs with a given tag
    */
@@ -93,8 +96,10 @@ class EigenProblem : public FEProblemBase
   unsigned int _n_eigen_pairs_required;
   bool _generalized_eigenvalue_problem;
   std::shared_ptr<NonlinearEigenSystem> _nl_eigen;
+
   bool _negative_sign_eigen_kernel;
 
+  unsigned int _active_eigen_index;
   /// Timers
   PerfID _compute_jacobian_tag_timer;
   PerfID _compute_jacobian_ab_timer;

diff --git a/framework/include/systems/NonlinearEigenSystem.h b/framework/include/systems/NonlinearEigenSystem.h
@@ -54,6 +54,16 @@ class NonlinearEigenSystem : public NonlinearSystemBase
 
   virtual NumericVector<Number> & RHS() override;
 
+  /*
+   *  A residual vector at MOOSE side for AX
+   */
+  NumericVector<Number> & ResidualVectorAX();
+
+  /*
+   * A residual vector at MOOSE side for BX
+   */
+  NumericVector<Number> & ResidualVectorBX();
+
   /**
    * Add the eigen tag to the right kernels
    */
@@ -141,6 +151,8 @@ class NonlinearEigenSystem : public NonlinearSystemBase
   EigenProblem & _eigen_problem;
   std::vector<std::pair<Real, Real>> _eigen_values;
   unsigned int _n_eigen_pairs_required;
+  NumericVector<Number> & _work_rhs_vector_AX;
+  NumericVector<Number> & _work_rhs_vector_BX;
   TagID _Ax_tag;
   TagID _Bx_tag;
   TagID _A_tag;
@@ -164,4 +176,3 @@ class NonlinearEigenSystem : public libMesh::ParallelObject
 };
 
 #endif
-
diff --git a/framework/src/executioners/Eigenvalue.C b/framework/src/executioners/Eigenvalue.C
@@ -43,7 +43,7 @@ Eigenvalue::Eigenvalue(const InputParameters & parameters)
 {
 // Extract and store SLEPc options
 #if LIBMESH_HAVE_SLEPC
-  Moose::SlepcSupport::storeSlepcOptions(_fe_problem, _pars);
+  Moose::SlepcSupport::storeSlepcOptions(_fe_problem, parameters);
 
   Moose::SlepcSupport::storeSlepcEigenProblemOptions(_eigen_problem, parameters);
   _eigen_problem.setEigenproblemType(_eigen_problem.solverParams()._eigen_problem_type);
@@ -66,5 +66,8 @@ Eigenvalue::execute()
   }
 #endif
 #endif
+
+  std::cout<<"Eigenvalue::execute()"<<std::endl;
+
   Steady::execute();
 }
diff --git a/framework/src/problems/EigenProblem.C b/framework/src/problems/EigenProblem.C
@@ -36,6 +36,9 @@ validParams<EigenProblem>()
                         "Whether or not to use a negative sign for eigenvalue kernels. "
                         "Using a negative sign makes eigenvalue kernels consistent with "
                         "a nonlinear solver");
+
+  params.addParam<unsigned int>("active_eigen_index", 0, "Which eigen vector is used to compute residual and also associateed to nonlinear variable");
+
   return params;
 }
 
@@ -46,6 +49,7 @@ EigenProblem::EigenProblem(const InputParameters & parameters)
     _generalized_eigenvalue_problem(false),
     _nl_eigen(std::make_shared<NonlinearEigenSystem>(*this, "eigen0")),
     _negative_sign_eigen_kernel(getParam<bool>("negative_sign_eigen_kernel")),
+    _active_eigen_index(getParam<unsigned int>("active_eigen_index")),
     _compute_jacobian_tag_timer(registerTimedSection("computeJacobianTag", 3)),
     _compute_jacobian_ab_timer(registerTimedSection("computeJacobianAB", 3)),
     _compute_residual_tag_timer(registerTimedSection("computeResidualTag", 3)),
@@ -61,7 +65,6 @@ EigenProblem::EigenProblem(const InputParameters & parameters)
   FEProblemBase::initNullSpaceVectors(parameters, *_nl_eigen);
 
   _eq.parameters.set<EigenProblem *>("_eigen_problem") = this;
-
 #else
   mooseError("Need to install SLEPc to solve eigenvalue problems, please reconfigure\n");
 #endif /* LIBMESH_HAVE_SLEPC */
@@ -214,6 +217,25 @@ EigenProblem::computeResidualAB(const NumericVector<Number> & soln,
   _nl_eigen->disassociateVectorFromTag(residualB, tagB);
 }
 
+Real
+EigenProblem::computeResidualL2Norm()
+{
+  computeResidualAB(*_nl_eigen->currentSolution(), _nl_eigen->ResidualVectorAX(), _nl_eigen->ResidualVectorBX(), _nl_eigen->nonEigenVectorTag(), _nl_eigen->eigenVectorTag());
+
+  Real eigenvalue = 1.0;
+
+  if (_active_eigen_index<_nl_eigen->getNumConvergedEigenvalues())
+    eigenvalue = _nl_eigen->getNthConvergedEigenvalue(_active_eigen_index).first;
+
+  // Scale BX with eigenvalue
+  _nl_eigen->ResidualVectorBX() *= eigenvalue;
+
+  // Compute entire residual
+  _nl_eigen->ResidualVectorAX() -= _nl_eigen->ResidualVectorBX();
+
+  return _nl_eigen->ResidualVectorAX().l2_norm();
+}
+
 #endif
 
 void

diff --git a/framework/src/systems/NonlinearEigenSystem.C b/framework/src/systems/NonlinearEigenSystem.C
@@ -129,7 +129,9 @@ NonlinearEigenSystem::NonlinearEigenSystem(EigenProblem & eigen_problem, const s
         eigen_problem, eigen_problem.es().add_system<TransientEigenSystem>(name), name),
     _transient_sys(eigen_problem.es().get_system<TransientEigenSystem>(name)),
     _eigen_problem(eigen_problem),
-    _n_eigen_pairs_required(eigen_problem.getNEigenPairsRequired())
+    _n_eigen_pairs_required(eigen_problem.getNEigenPairsRequired()),
+    _work_rhs_vector_AX(addVector("work_rhs_vector_Ax",false,PARALLEL)),
+    _work_rhs_vector_BX(addVector("work_rhs_vector_Bx",false,PARALLEL))
 {
   sys().attach_assemble_function(Moose::assemble_matrix);
 
@@ -248,8 +250,19 @@ NonlinearEigenSystem::getCurrentNonlinearIterationNumber()
 NumericVector<Number> &
 NonlinearEigenSystem::RHS()
 {
-  mooseError("did not implement yet \n");
-  // return NULL;
+  return _work_rhs_vector_BX;
+}
+
+NumericVector<Number> &
+NonlinearEigenSystem::ResidualVectorAX()
+{
+  return _work_rhs_vector_AX;
+}
+
+NumericVector<Number> &
+NonlinearEigenSystem::ResidualVectorBX()
+{
+  return _work_rhs_vector_BX;
 }
 
 NonlinearSolver<Number> *

diff --git a/test/tests/kernels/2d_diffusion/2d_diffusion_system.i b/test/tests/kernels/2d_diffusion/2d_diffusion_system.i
@@ -0,0 +1,78 @@
+###########################################################
+# This is a simple test of the Kernel System.
+# It solves the Laplacian equation on a small 2x2 grid.
+# The "Diffusion" kernel is used to calculate the
+# residuals of the weak form of this operator.
+#
+# @Requirement F3.30
+###########################################################
+
+[Mesh]
+  file = square.e
+[]
+
+[Variables]
+  [./u]
+    order = FIRST
+    family = LAGRANGE
+  [../]
+
+  [./v]
+    order = FIRST
+    family = LAGRANGE
+  [../]
+[]
+
+[Kernels]
+  [./diff]
+    type = Diffusion
+    variable = u
+  [../]
+  [./diff1]
+    type = Diffusion
+    variable = v
+  [../]
+[]
+
+[BCs]
+  active = 'left right'
+
+  [./left]
+    type = DirichletBC
+    variable = u
+    boundary = 1
+    value = 0
+  [../]
+
+  [./right]
+    type = DirichletBC
+    variable = u
+    boundary = 2
+    value = 1
+  [../]
+
+  [./left1]
+    type = DirichletBC
+    variable = v
+    boundary = 1
+    value = 0
+  [../]
+
+  [./right1]
+    type = DirichletBC
+    variable = v
+    boundary = 2
+    value = 10
+  [../]
+[]
+
+[Executioner]
+  type = Steady
+
+  solve_type = 'NEWTON'
+[]
+
+[Outputs]
+  file_base = out
+  exodus = true
+[]
diff --git a/test/tests/preconditioners/hyprematrix/2d_diffusion_hyprematrix.i b/test/tests/preconditioners/hyprematrix/2d_diffusion_hyprematrix.i
@@ -0,0 +1,56 @@
+[Mesh]
+  type = GeneratedMesh
+  dim = 2
+  nx = 4
+  ny = 4
+[]
+
+[Variables]
+  active = 'u'
+
+  [./u]
+    order = FIRST
+    family = LAGRANGE
+  [../]
+[]
+
+[Kernels]
+  active = 'diff'
+
+  [./diff]
+    type = Diffusion
+    variable = u
+  [../]
+[]
+
+[BCs]
+  active = 'left right'
+
+  [./left]
+    type = DirichletBC
+    variable = u
+    boundary = 1
+    value = 0
+  [../]
+
+  [./right]
+    type = DirichletBC
+    variable = u
+    boundary = 2
+    value = 1
+  [../]
+[]
+
+[Executioner]
+  type = Steady
+  solve_type = 'NEWTON'
+
+  petsc_options_iname = '-pc_type'
+  petsc_options_value = 'hypre'
+
+  hypre_matrix = true
+[]
+
+[Outputs]
+  exodus = true
+[]
diff --git a/test/tests/problems/eigen_problem/ne_coupled_picard_masterT.i b/test/tests/problems/eigen_problem/ne_coupled_picard_masterT.i
@@ -49,7 +49,7 @@
 [Executioner]
   type = Steady
   picard_max_its = 30
-  picard_rel_tol = 1e-60
+  picard_rel_tol = 1e-6
 []
 
 [MultiApps]

diff --git a/test/tests/problems/eigen_problem/ne_coupled_picard_masterT_sub.i b/test/tests/problems/eigen_problem/ne_coupled_picard_masterT_sub.i
@@ -49,7 +49,6 @@
     variable = power
     source_variable = u
     normalization = unorm
-    # this coefficient will affect the eigenvalue.
     normal_factor = 10
     execute_on = timestep_end
   [../]

diff --git a/test/tests/problems/eigen_problem/tests b/test/tests/problems/eigen_problem/tests
@@ -111,4 +111,6 @@
     slepc_version = '>=3.8.0'
     requirement = "Eigenvalue system should be able to handle DG kernels"
   []
+
+
 []