www.mooseframework.org
NonlinearEigenSystem.C
Go to the documentation of this file.
1 /****************************************************************/
2 /* DO NOT MODIFY THIS HEADER */
3 /* MOOSE - Multiphysics Object Oriented Simulation Environment */
4 /* */
5 /* (c) 2010 Battelle Energy Alliance, LLC */
6 /* ALL RIGHTS RESERVED */
7 /* */
8 /* Prepared by Battelle Energy Alliance, LLC */
9 /* Under Contract No. DE-AC07-05ID14517 */
10 /* With the U. S. Department of Energy */
11 /* */
12 /* See COPYRIGHT for full restrictions */
13 /****************************************************************/
14 
15 #include "NonlinearEigenSystem.h"
16 
17 // MOOSE includes
18 #include "DirichletBC.h"
19 #include "EigenProblem.h"
20 #include "IntegratedBC.h"
21 #include "KernelBase.h"
22 #include "NodalBC.h"
23 #include "TimeIntegrator.h"
24 #include "SlepcSupport.h"
25 
26 #include "libmesh/eigen_system.h"
27 #include "libmesh/libmesh_config.h"
28 #include "libmesh/petsc_matrix.h"
29 #include "libmesh/sparse_matrix.h"
30 
31 #if LIBMESH_HAVE_SLEPC
32 
33 namespace Moose
34 {
35 
36 void
37 assemble_matrix(EquationSystems & es, const std::string & system_name)
38 {
39  EigenProblem * p = es.parameters.get<EigenProblem *>("_eigen_problem");
40  EigenSystem & eigen_system = es.get_system<EigenSystem>(system_name);
41 
43  {
44  p->computeJacobian(
45  *eigen_system.current_local_solution.get(), *eigen_system.matrix_A, Moose::KT_NONEIGEN);
46  }
47  else
48  {
49  Mat petsc_mat_A = static_cast<PetscMatrix<Number> &>(*eigen_system.matrix_A).mat();
50 
51  PetscObjectComposeFunction((PetscObject)petsc_mat_A,
52  "formJacobian",
54  PetscObjectComposeFunction((PetscObject)petsc_mat_A,
55  "formFunction",
57 
58  PetscContainer container;
59  PetscContainerCreate(eigen_system.comm().get(), &container);
60  PetscContainerSetPointer(container, p);
61  PetscObjectCompose((PetscObject)petsc_mat_A, "formJacobianCtx", nullptr);
62  PetscObjectCompose((PetscObject)petsc_mat_A, "formJacobianCtx", (PetscObject)container);
63  PetscObjectCompose((PetscObject)petsc_mat_A, "formFunctionCtx", nullptr);
64  PetscObjectCompose((PetscObject)petsc_mat_A, "formFunctionCtx", (PetscObject)container);
65  PetscContainerDestroy(&container);
66  }
67  if (eigen_system.generalized())
68  {
69  if (eigen_system.matrix_B)
70  {
72  {
73  p->computeJacobian(
74  *eigen_system.current_local_solution.get(), *eigen_system.matrix_B, Moose::KT_EIGEN);
75  }
76  else
77  {
78  Mat petsc_mat_B = static_cast<PetscMatrix<Number> &>(*eigen_system.matrix_B).mat();
79 
80  PetscObjectComposeFunction((PetscObject)petsc_mat_B,
81  "formJacobian",
83  PetscObjectComposeFunction((PetscObject)petsc_mat_B,
84  "formFunction",
86 
87  PetscContainer container;
88  PetscContainerCreate(eigen_system.comm().get(), &container);
89  PetscContainerSetPointer(container, p);
90  PetscObjectCompose((PetscObject)petsc_mat_B, "formFunctionCtx", nullptr);
91  PetscObjectCompose((PetscObject)petsc_mat_B, "formFunctionCtx", (PetscObject)container);
92  PetscObjectCompose((PetscObject)petsc_mat_B, "formJacobianCtx", nullptr);
93  PetscObjectCompose((PetscObject)petsc_mat_B, "formJacobianCtx", (PetscObject)container);
94  PetscContainerDestroy(&container);
95  }
96  }
97  else
98  mooseError("It is a generalized eigenvalue problem but matrix B is empty\n");
99  }
100 }
101 }
102 
103 NonlinearEigenSystem::NonlinearEigenSystem(EigenProblem & eigen_problem, const std::string & name)
105  eigen_problem, eigen_problem.es().add_system<TransientEigenSystem>(name), name),
106  _transient_sys(eigen_problem.es().get_system<TransientEigenSystem>(name)),
107  _eigen_problem(eigen_problem),
108  _n_eigen_pairs_required(eigen_problem.getNEigenPairsRequired())
109 {
110  sys().attach_assemble_function(Moose::assemble_matrix);
111 }
112 
113 void
115 {
116  // Clear the iteration counters
117  _current_l_its.clear();
118  _current_nl_its = 0;
119  // Initialize the solution vector using a predictor and known values from nodal bcs
121 
122  // Solve the transient problem if we have a time integrator; the
123  // steady problem if not.
124  if (_time_integrator)
125  {
126  _time_integrator->solve();
127  _time_integrator->postSolve();
128  }
129  else
130  system().solve();
131 
132  // store eigenvalues
133  unsigned int n_converged_eigenvalues = getNumConvergedEigenvalues();
134 
136 
137  if (_n_eigen_pairs_required < n_converged_eigenvalues)
138  n_converged_eigenvalues = _n_eigen_pairs_required;
139 
140  _eigen_values.resize(n_converged_eigenvalues);
141  for (unsigned int n = 0; n < n_converged_eigenvalues; n++)
143 }
144 
145 void
147 {
148  mooseError("did not implement yet \n");
149 }
150 
151 void
153 {
154  mooseError("did not implement yet \n");
155 }
156 
157 bool
159 {
160  return _transient_sys.get_n_converged();
161 }
162 
163 unsigned int
165 {
166  mooseError("did not implement yet \n");
167  return 0;
168 }
169 
170 NumericVector<Number> &
172 {
173  mooseError("did not implement yet \n");
174  // return NULL;
175 }
176 
177 NonlinearSolver<Number> *
179 {
180  mooseError("did not implement yet \n");
181  return NULL;
182 }
183 
184 void
185 NonlinearEigenSystem::addEigenKernels(std::shared_ptr<KernelBase> kernel, THREAD_ID tid)
186 {
187  if (kernel->isEigenKernel())
188  _eigen_kernels.addObject(kernel, tid);
189  else
190  _non_eigen_kernels.addObject(kernel, tid);
191 }
192 
193 void
195 {
197  mooseError("Can't set an inhomogeneous integrated boundary condition for eigenvalue problems.");
198 
200  {
201  const auto & nodal_bcs = _nodal_bcs.getActiveObjects();
202  for (const auto & nodal_bc : nodal_bcs)
203  {
204  std::shared_ptr<DirichletBC> nbc = std::dynamic_pointer_cast<DirichletBC>(nodal_bc);
205  if (nbc && nbc->getParam<Real>("value"))
206  mooseError(
207  "Can't set an inhomogeneous Dirichlet boundary condition for eigenvalue problems.");
208  else if (!nbc)
209  mooseError("Invalid NodalBC for eigenvalue problems, please use homogeneous Dirichlet.");
210  }
211  }
212 }
213 
214 const std::pair<Real, Real>
216 {
217  unsigned int n_converged_eigenvalues = getNumConvergedEigenvalues();
218  if (n >= n_converged_eigenvalues)
219  mooseError(n, " not in [0, ", n_converged_eigenvalues, ")");
220  return _transient_sys.get_eigenpair(n);
221 }
222 
223 #else
224 
226  const std::string & /*name*/)
227  : libMesh::ParallelObject(eigen_problem)
228 {
229  mooseError("Need to install SLEPc to solve eigenvalue problems, please reconfigure libMesh\n");
230 }
231 
232 #endif /* LIBMESH_HAVE_SLEPC */
virtual unsigned int getNEigenPairsRequired()
Definition: EigenProblem.h:43
virtual bool isNonlinearEigenvalueSolver()
Definition: EigenProblem.C:155
virtual TransientEigenSystem & sys()
virtual void addEigenKernels(std::shared_ptr< KernelBase > kernel, THREAD_ID tid) override
virtual NonlinearSolver< Number > * nonlinearSolver() override
PetscErrorCode mooseSlepcEigenFormJacobianB(SNES snes, Vec x, Mat jac, Mat pc, void *ctx)
Definition: SlepcSupport.C:488
void mooseError(Args &&...args)
Emit an error message with the given stringified, concatenated args and terminate the application...
Definition: MooseError.h:182
MooseObjectWarehouse< NodalBC > _nodal_bcs
KernelWarehouse _eigen_kernels
PetscErrorCode mooseSlepcEigenFormFunctionA(SNES snes, Vec x, Vec r, void *ctx)
Definition: SlepcSupport.C:516
std::shared_ptr< TimeIntegrator > _time_integrator
Time integrator.
NonlinearEigenSystem(EigenProblem &problem, const std::string &name)
Nonlinear system to be solved.
std::vector< unsigned int > _current_l_its
virtual const std::pair< Real, Real > getNthConvergedEigenvalue(dof_id_type n)
Return the Nth converged eigenvalue.
virtual unsigned int getNumConvergedEigenvalues() const
Get the number of converged eigenvalues.
virtual bool converged() override
Returns the convergence state.
Boundary condition of a Dirichlet type.
Definition: DirichletBC.h:30
KernelWarehouse _non_eigen_kernels
bool hasActiveObjects(THREAD_ID tid=0) const
Convenience functions for determining if objects exist.
virtual void stopSolve() override
Quit the current solve as soon as possible.
void assemble_matrix(EquationSystems &es, const std::string &system_name)
PetscErrorCode mooseSlepcEigenFormJacobianA(SNES snes, Vec x, Mat jac, Mat pc, void *ctx)
Definition: SlepcSupport.C:479
virtual void setupFiniteDifferencedPreconditioner() override
PetscErrorCode mooseSlepcEigenFormFunctionB(SNES snes, Vec x, Vec r, void *ctx)
Definition: SlepcSupport.C:525
virtual unsigned int getCurrentNonlinearIterationNumber() override
Returns the current nonlinear iteration number.
TransientEigenSystem & _transient_sys
virtual void computeJacobian(const NumericVector< Number > &soln, SparseMatrix< Number > &jacobian, Moose::KernelType kernel_type) override
Definition: EigenProblem.C:115
virtual System & system() override
Get the reference to the libMesh system.
PetscInt n
void addObject(std::shared_ptr< KernelBase > object, THREAD_ID tid=0) override
Add Kernel to the storage structure.
const InputParameters & parameters() const
Get the parameters of the object.
Definition: MooseObject.h:53
MooseObjectWarehouse< IntegratedBC > _integrated_bcs
Definition: Moose.h:84
Specialization of SubProblem for solving nonlinear equations plus auxiliary equations.
Definition: EigenProblem.h:32
std::vector< std::pair< Real, Real > > _eigen_values
virtual NumericVector< Number > & RHS() override
unsigned int _n_eigen_pairs_required
virtual void solve() override
Solve the system (using libMesh magic)
const std::vector< std::shared_ptr< T > > & getActiveObjects(THREAD_ID tid=0) const
Retrieve complete vector to the active all/block/boundary restricted objects for a given thread...
void checkIntegrity()
For eigenvalue problems (including standard and generalized), inhomogeneous (Dirichlet or Neumann) bo...
unsigned int THREAD_ID
Definition: MooseTypes.h:79