www.mooseframework.org
NonlinearEigenSystem.C
Go to the documentation of this file.
1 //* This file is part of the MOOSE framework
2 //* https://www.mooseframework.org
3 //*
4 //* All rights reserved, see COPYRIGHT for full restrictions
5 //* https://github.com/idaholab/moose/blob/master/COPYRIGHT
6 //*
7 //* Licensed under LGPL 2.1, please see LICENSE for details
8 //* https://www.gnu.org/licenses/lgpl-2.1.html
9 
10 #include "NonlinearEigenSystem.h"
11 
12 // MOOSE includes
13 #include "DirichletBC.h"
14 #include "EigenDirichletBC.h"
15 #include "EigenProblem.h"
16 #include "IntegratedBC.h"
17 #include "KernelBase.h"
18 #include "NodalBC.h"
19 #include "TimeIntegrator.h"
20 #include "SlepcSupport.h"
21 
22 #include "libmesh/eigen_system.h"
23 #include "libmesh/libmesh_config.h"
24 #include "libmesh/petsc_matrix.h"
25 #include "libmesh/sparse_matrix.h"
26 
27 #if LIBMESH_HAVE_SLEPC
28 
29 namespace Moose
30 {
31 
32 void
33 assemble_matrix(EquationSystems & es, const std::string & system_name)
34 {
35  EigenProblem * p = es.parameters.get<EigenProblem *>("_eigen_problem");
36  EigenSystem & eigen_system = es.get_system<EigenSystem>(system_name);
38 
39  // If it is a linear generalized eigenvalue problem,
40  // we assemble A and B together
41  if (!p->isNonlinearEigenvalueSolver() && eigen_system.generalized())
42  {
43  p->computeJacobianAB(*eigen_system.current_local_solution.get(),
44  *eigen_system.matrix_A,
45  *eigen_system.matrix_B,
46  eigen_nl.nonEigenMatrixTag(),
47  eigen_nl.eigenMatrixTag());
48  return;
49  }
50 
52  {
53  p->computeJacobianTag(*eigen_system.current_local_solution.get(),
54  *eigen_system.matrix_A,
55  eigen_nl.nonEigenMatrixTag());
56  }
57  else
58  {
59  Mat petsc_mat_A = static_cast<PetscMatrix<Number> &>(*eigen_system.matrix_A).mat();
60 
61  PetscObjectComposeFunction((PetscObject)petsc_mat_A,
62  "formJacobian",
64  PetscObjectComposeFunction((PetscObject)petsc_mat_A,
65  "formFunction",
67 
68  PetscObjectComposeFunction((PetscObject)petsc_mat_A,
69  "formFunctionAB",
71 
72  PetscContainer container;
73  PetscContainerCreate(eigen_system.comm().get(), &container);
74  PetscContainerSetPointer(container, p);
75  PetscObjectCompose((PetscObject)petsc_mat_A, "formJacobianCtx", nullptr);
76  PetscObjectCompose((PetscObject)petsc_mat_A, "formJacobianCtx", (PetscObject)container);
77  PetscObjectCompose((PetscObject)petsc_mat_A, "formFunctionCtx", nullptr);
78  PetscObjectCompose((PetscObject)petsc_mat_A, "formFunctionCtx", (PetscObject)container);
79  PetscContainerDestroy(&container);
80 
81  // Let libmesh do not close matrices before solve
82  eigen_system.eigen_solver->set_close_matrix_before_solve(false);
83  }
84  if (eigen_system.generalized())
85  {
86  if (eigen_system.matrix_B)
87  {
89  {
90  p->computeJacobianTag(*eigen_system.current_local_solution.get(),
91  *eigen_system.matrix_B,
92  eigen_nl.eigenMatrixTag());
93  }
94  else
95  {
96  Mat petsc_mat_B = static_cast<PetscMatrix<Number> &>(*eigen_system.matrix_B).mat();
97 
98  PetscObjectComposeFunction((PetscObject)petsc_mat_B,
99  "formJacobian",
101  PetscObjectComposeFunction((PetscObject)petsc_mat_B,
102  "formFunction",
104 
105  PetscContainer container;
106  PetscContainerCreate(eigen_system.comm().get(), &container);
107  PetscContainerSetPointer(container, p);
108  PetscObjectCompose((PetscObject)petsc_mat_B, "formFunctionCtx", nullptr);
109  PetscObjectCompose((PetscObject)petsc_mat_B, "formFunctionCtx", (PetscObject)container);
110  PetscObjectCompose((PetscObject)petsc_mat_B, "formJacobianCtx", nullptr);
111  PetscObjectCompose((PetscObject)petsc_mat_B, "formJacobianCtx", (PetscObject)container);
112  PetscContainerDestroy(&container);
113  }
114  }
115  else
116  mooseError("It is a generalized eigenvalue problem but matrix B is empty\n");
117  }
118 }
119 }
120 
121 NonlinearEigenSystem::NonlinearEigenSystem(EigenProblem & eigen_problem, const std::string & name)
123  eigen_problem, eigen_problem.es().add_system<TransientEigenSystem>(name), name),
124  _transient_sys(eigen_problem.es().get_system<TransientEigenSystem>(name)),
125  _eigen_problem(eigen_problem),
126  _n_eigen_pairs_required(eigen_problem.getNEigenPairsRequired())
127 {
128  sys().attach_assemble_function(Moose::assemble_matrix);
129 
130  _Ax_tag = eigen_problem.addVectorTag("Ax_tag");
131 
132  _Bx_tag = eigen_problem.addVectorTag("Eigen");
133 
134  _A_tag = eigen_problem.addMatrixTag("A_tag");
135 
136  _B_tag = eigen_problem.addMatrixTag("Eigen");
137 }
138 
139 void
141 {
143 
145 
147 }
148 
149 template <typename T>
150 void
152 {
153  for (THREAD_ID tid = 0; tid < warehouse.numThreads(); tid++)
154  {
155  auto & objects = warehouse.getObjects(tid);
156 
157  for (auto & object : objects)
158  {
159  auto & vtags = object->getVectorTags();
160  // If this is not an eigen kernel
161  if (vtags.find(_Bx_tag) == vtags.end())
162  object->useVectorTag(_Ax_tag);
163  else // also associate eigen matrix tag if this is a eigen kernel
164  object->useMatrixTag(_B_tag);
165 
166  auto & mtags = object->getMatrixTags();
167  if (mtags.find(_B_tag) == mtags.end())
168  object->useMatrixTag(_A_tag);
169  else
170  object->useVectorTag(_Bx_tag);
171  }
172  }
173 }
174 
175 void
177 {
178  // Clear the iteration counters
179  _current_l_its.clear();
180  _current_nl_its = 0;
181  // Initialize the solution vector using a predictor and known values from nodal bcs
183 
184 // In DEBUG mode, Libmesh will check the residual automatically. This may cause
185 // an error because B does not need to assembly by default.
186 #ifdef DEBUG
188  sys().matrix_B->close();
189 #endif
190  // Solve the transient problem if we have a time integrator; the
191  // steady problem if not.
192  if (_time_integrator)
193  {
194  _time_integrator->solve();
195  _time_integrator->postSolve();
196  }
197  else
198  system().solve();
199 
200  // store eigenvalues
201  unsigned int n_converged_eigenvalues = getNumConvergedEigenvalues();
202 
204 
205  if (_n_eigen_pairs_required < n_converged_eigenvalues)
206  n_converged_eigenvalues = _n_eigen_pairs_required;
207 
208  _eigen_values.resize(n_converged_eigenvalues);
209  for (unsigned int n = 0; n < n_converged_eigenvalues; n++)
211 }
212 
213 void
215 {
216  mooseError("did not implement yet \n");
217 }
218 
219 void
221 {
222  mooseError("did not implement yet \n");
223 }
224 
225 bool
227 {
228  return _transient_sys.get_n_converged();
229 }
230 
231 unsigned int
233 {
234  mooseError("did not implement yet \n");
235  return 0;
236 }
237 
238 NumericVector<Number> &
240 {
241  mooseError("did not implement yet \n");
242  // return NULL;
243 }
244 
245 NonlinearSolver<Number> *
247 {
248  mooseError("did not implement yet \n");
249  return NULL;
250 }
251 
252 void
254 {
256  mooseError("Can't set an inhomogeneous integrated boundary condition for eigenvalue problems.");
257 
259  {
260  const auto & nodal_bcs = _nodal_bcs.getActiveObjects();
261  for (const auto & nodal_bc : nodal_bcs)
262  {
263  auto nbc = std::dynamic_pointer_cast<DirichletBC>(nodal_bc);
264  auto eigen_nbc = std::dynamic_pointer_cast<EigenDirichletBC>(nodal_bc);
265  if (nbc && nbc->getParam<Real>("value"))
266  mooseError(
267  "Can't set an inhomogeneous Dirichlet boundary condition for eigenvalue problems.");
268  else if (!nbc && !eigen_nbc)
269  mooseError("Invalid NodalBC for eigenvalue problems, please use homogeneous Dirichlet.");
270  }
271  }
272 }
273 
274 const std::pair<Real, Real>
276 {
277  unsigned int n_converged_eigenvalues = getNumConvergedEigenvalues();
278  if (n >= n_converged_eigenvalues)
279  mooseError(n, " not in [0, ", n_converged_eigenvalues, ")");
280  return _transient_sys.get_eigenpair(n);
281 }
282 
283 #else
284 
286  const std::string & /*name*/)
287  : libMesh::ParallelObject(eigen_problem)
288 {
289  mooseError("Need to install SLEPc to solve eigenvalue problems, please reconfigure libMesh\n");
290 }
291 
292 #endif /* LIBMESH_HAVE_SLEPC */
virtual unsigned int getNEigenPairsRequired()
Definition: EigenProblem.h:34
Nonlinear eigenvalue system to be solved.
virtual bool isNonlinearEigenvalueSolver()
Definition: EigenProblem.C:243
virtual TagID addVectorTag(TagName tag_name)
Create a Tag.
Definition: SubProblem.C:66
virtual TransientEigenSystem & sys()
THREAD_ID numThreads()
Return the number of threads.
void mooseError(Args &&... args)
Emit an error message with the given stringified, concatenated args and terminate the application...
Definition: MooseError.h:207
virtual NonlinearSolver< Number > * nonlinearSolver() override
PetscErrorCode mooseSlepcEigenFormJacobianB(SNES snes, Vec x, Mat jac, Mat pc, void *ctx)
Definition: SlepcSupport.C:480
void addEigenTagToMooseObjects(MooseObjectTagWarehouse< T > &warehouse)
Add the eigen tag to the right kernels.
MooseObjectTagWarehouse< NodalBCBase > _nodal_bcs
PetscErrorCode mooseSlepcEigenFormFunctionAB(SNES snes, Vec x, Vec Ax, Vec Bx, void *ctx)
Definition: SlepcSupport.C:535
std::unique_ptr< T_DEST, T_DELETER > dynamic_pointer_cast(std::unique_ptr< T_SRC, T_DELETER > &src)
These are reworked from https://stackoverflow.com/a/11003103.
virtual void initialSetup() override
PetscErrorCode mooseSlepcEigenFormFunctionA(SNES snes, Vec x, Vec r, void *ctx)
Definition: SlepcSupport.C:511
The following methods are specializations for using the libMesh::Parallel::packed_range_* routines fo...
virtual unsigned int getNumConvergedEigenvalues() const
Get the number of converged eigenvalues.
NonlinearEigenSystem(EigenProblem &problem, const std::string &name)
Nonlinear system to be solved.
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...
std::vector< unsigned int > _current_l_its
const std::vector< std::shared_ptr< T > > & getObjects(THREAD_ID tid=0) const
Retrieve complete vector to the all/block/boundary restricted objects for a given thread...
virtual const std::pair< Real, Real > getNthConvergedEigenvalue(dof_id_type n)
Return the Nth converged eigenvalue.
const InputParameters & parameters() const
Get the parameters of the object.
Definition: MooseObject.h:65
virtual bool converged() override
Returns the convergence state.
virtual TagID addMatrixTag(TagName tag_name)
Create a Tag.
Definition: SubProblem.C:112
Boundary condition of a Dirichlet type.
Definition: DirichletBC.h:24
MooseObjectTagWarehouse< KernelBase > _kernels
virtual void initialSetup()
virtual void stopSolve() override
Quit the current solve as soon as possible.
void assemble_matrix(EquationSystems &es, const std::string &system_name)
A storage container for MooseObjects that inherit from SetupInterface.
PetscErrorCode mooseSlepcEigenFormJacobianA(SNES snes, Vec x, Mat jac, Mat pc, void *ctx)
Definition: SlepcSupport.C:468
Set Dirichlet boundary condition for eigenvalue problems.
virtual void setupFiniteDifferencedPreconditioner() override
PetscErrorCode mooseSlepcEigenFormFunctionB(SNES snes, Vec x, Vec r, void *ctx)
Definition: SlepcSupport.C:523
virtual unsigned int getCurrentNonlinearIterationNumber() override
Returns the current nonlinear iteration number.
TransientEigenSystem & _transient_sys
bool hasActiveObjects(THREAD_ID tid=0) const
virtual System & system() override
Get the reference to the libMesh system.
PetscInt n
std::shared_ptr< TimeIntegrator > _time_integrator
Time integrator.
Definition: SystemBase.h:782
Definition: Moose.h:112
Problem for solving eigenvalue problems.
Definition: EigenProblem.h:25
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)
MooseObjectTagWarehouse< IntegratedBCBase > _integrated_bcs
virtual bool isGeneralizedEigenvalueProblem()
Definition: EigenProblem.h:39
void checkIntegrity()
For eigenvalue problems (including standard and generalized), inhomogeneous (Dirichlet or Neumann) bo...
NonlinearEigenSystem & getNonlinearEigenSystem()
Definition: EigenProblem.h:44
virtual void computeJacobianAB(const NumericVector< Number > &soln, SparseMatrix< Number > &jacobianA, SparseMatrix< Number > &jacobianB, TagID tagA, TagID tagB)
Form two Jacobian matrices, whre each is associateed with one tag, through one element-loop.
Definition: EigenProblem.C:134
virtual void computeJacobianTag(const NumericVector< Number > &soln, SparseMatrix< Number > &jacobian, TagID tag) override
Form a Jacobian matrix for all kernels and BCs with a given tag.
Definition: EigenProblem.C:112
unsigned int THREAD_ID
Definition: MooseTypes.h:161