www.mooseframework.org
EigenExecutionerBase.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 "EigenExecutionerBase.h"
16 
17 // MOOSE includes
18 #include "AuxiliarySystem.h"
19 #include "DisplacedProblem.h"
20 #include "FEProblem.h"
21 #include "MooseApp.h"
22 #include "MooseEigenSystem.h"
23 #include "UserObject.h"
24 
25 template <>
28 {
30  params.addRequiredParam<PostprocessorName>("bx_norm", "To evaluate |Bx| for the eigenvalue");
31  params.addParam<PostprocessorName>("normalization", "To evaluate |x| for normalization");
32  params.addParam<Real>("normal_factor", "Normalize x to make |x| equal to this factor");
33  params.addParam<bool>(
34  "output_before_normalization", true, "True to output a step before normalization");
35  params.addParam<bool>("auto_initialization", true, "True to ask the solver to set initial");
36  params.addParam<Real>("time", 0.0, "System time");
37 
38  params.addPrivateParam<bool>("_eigen", true);
39 
40  params.addParamNamesToGroup("normalization normal_factor output_before_normalization",
41  "Normalization");
42  params.addParamNamesToGroup("auto_initialization time", "Advanced");
43 
44  params.addPrivateParam<bool>("_eigen", true);
45 
46  return params;
47 }
48 
49 const Real &
51 {
52  return _source_integral_old;
53 }
54 
56  : Executioner(parameters),
58  _eigen_sys(static_cast<MooseEigenSystem &>(_problem.getNonlinearSystemBase())),
59  _eigenvalue(declareRestartableData("eigenvalue", 1.0)),
62  _normalization(isParamValid("normalization")
63  ? getPostprocessorValue("normalization")
64  : getPostprocessorValue("bx_norm")) // use |Bx| for normalization by default
65 {
66  // FIXME: currently we have to use old and older solution vectors for power iteration.
67  // We will need 'step' in the future.
68  _problem.transient(true);
69 
70  // we want to tell the App about what our system time is (in case anyone else is interested).
71  Real system_time = getParam<Real>("time");
72  _app.setStartTime(system_time);
73 
74  // set the system time
75  _problem.time() = system_time;
76  _problem.timeOld() = system_time;
77 
78  // used for controlling screen print-out
79  _problem.timeStep() = 0;
80  _problem.dt() = 1.0;
81 }
82 
83 void
85 {
88 
89  if (getParam<bool>("auto_initialization"))
90  {
91  // Initialize the solution of the eigen variables
92  // Note: initial conditions will override this if there is any by _problem.initialSetup()
94  }
97 
98  // check when the postprocessors are evaluated
99  ExecFlagType bx_execflag =
100  _problem.getUserObject<UserObject>(getParam<PostprocessorName>("bx_norm")).execBitFlags();
101  if ((bx_execflag & EXEC_LINEAR) == EXEC_NONE)
102  mooseError("Postprocessor " + getParam<PostprocessorName>("bx_norm") +
103  " requires execute_on = 'linear'");
104 
105  if (isParamValid("normalization"))
107  _problem.getUserObject<UserObject>(getParam<PostprocessorName>("normalization"))
108  .execBitFlags();
109  else
110  _norm_execflag = bx_execflag;
111 
112  // check if _source_integral has been evaluated during initialSetup()
113  if ((bx_execflag & EXEC_INITIAL) == EXEC_NONE)
114  _problem.execute(EXEC_LINEAR);
115 
116  if (_source_integral == 0.0)
117  mooseError("|Bx| = 0!");
118 
119  // normalize solution to make |Bx|=_eigenvalue, _eigenvalue at this point has the initialized
120  // value
122 
123  if (_problem.getDisplacedProblem() != NULL)
124  _problem.getDisplacedProblem()->syncSolutions();
125 
126  /* a time step check point */
128 }
129 
130 void
132 {
133  Real consistency_tolerance = 1e-10;
134 
135  // Scale the solution so that the postprocessor is equal to k.
136  // Note: all dependent objects of k must be evaluated on linear!
137  // We have a fix point loop here, in case the postprocessor is a nonlinear function of the scaling
138  // factor.
139  // FIXME: We have assumed this loop always converges.
140  while (std::fabs(k - _source_integral) > consistency_tolerance * std::fabs(k))
141  {
142  // On the first time entering, the _source_integral has been updated properly in
143  // FEProblemBase::initialSetup()
146  std::stringstream ss;
147  ss << std::fixed << std::setprecision(10) << _source_integral;
148  _console << "\n|Bx| = " << ss.str() << std::endl;
149  }
150 }
151 
152 void
154 {
155  // check to make sure that we don't have any time kernels in this simulation
157  mooseError("You have specified time kernels in your steady state eigenvalue simulation");
159  mooseError("You have not specified any eigen kernels in your eigenvalue simulation");
160 }
161 
162 void
164  unsigned int max_iter,
165  Real pfactor,
166  bool cheb_on,
167  Real tol_eig,
168  bool echo,
169  PostprocessorName xdiff,
170  Real tol_x,
171  Real & k,
172  Real & initial_res)
173 {
174  mooseAssert(max_iter >= min_iter,
175  "Maximum number of power iterations must be greater than or equal to its minimum");
176  mooseAssert(pfactor > 0.0, "Invaid linear convergence tolerance");
177  mooseAssert(tol_eig > 0.0, "Invalid eigenvalue tolerance");
178  mooseAssert(tol_x > 0.0, "Invalid solution norm tolerance");
179 
180  // obtain the solution diff
181  const PostprocessorValue * solution_diff = NULL;
182  if (xdiff != "")
183  {
184  solution_diff = &getPostprocessorValueByName(xdiff);
185  ExecFlagType xdiff_execflag = _problem.getUserObject<UserObject>(xdiff).execBitFlags();
186  if ((xdiff_execflag & EXEC_LINEAR) == EXEC_NONE)
187  mooseError("Postprocessor " + xdiff + " requires execute_on = 'linear'");
188  }
189 
190  // not perform any iteration when max_iter==0
191  if (max_iter == 0)
192  return;
193 
194  // turn off nonlinear flag so that RHS kernels opterate on previous solutions
196 
197  // FIXME: currently power iteration use old and older solutions,
198  // so save old and older solutions before they are changed by the power iteration
200  if (_problem.getDisplacedProblem() != NULL)
201  _problem.getDisplacedProblem()->saveOldSolutions();
202 
203  // save solver control parameters to be modified by the power iteration
204  Real tol1 = _problem.es().parameters.get<Real>("linear solver tolerance");
205  unsigned int num1 =
206  _problem.es().parameters.get<unsigned int>("nonlinear solver maximum iterations");
207  Real tol2 = _problem.es().parameters.get<Real>("nonlinear solver relative residual tolerance");
208 
209  // every power iteration is a linear solve, so set nonlinear iteration number to one
210  _problem.es().parameters.set<Real>("linear solver tolerance") = pfactor;
211  // disable nonlinear convergence check
212  _problem.es().parameters.set<unsigned int>("nonlinear solver maximum iterations") = 1;
213  _problem.es().parameters.set<Real>("nonlinear solver relative residual tolerance") = 1 - 1e-8;
214 
215  if (echo)
216  {
217  _console << '\n';
218  _console << " Power iterations starts\n";
219  _console << " ________________________________________________________________________________ "
220  << std::endl;
221  }
222 
223  // some iteration variables
224  Chebyshev_Parameters chebyshev_parameters;
225 
226  std::vector<Real> keff_history;
227  std::vector<Real> diff_history;
228 
229  unsigned int iter = 0;
230 
231  // power iteration loop...
232  // Note: |Bx|/k will stay constant one!
233  makeBXConsistent(k);
234  while (true)
235  {
236  if (echo)
237  _console << " Power iteration= " << iter << std::endl;
238 
239  // Important: we do not call _problem.advanceState() because we do not
240  // want to overwrite the old postprocessor values and old material
241  // properties in stateful materials.
244  if (_problem.getDisplacedProblem() != NULL)
245  {
246  _problem.getDisplacedProblem()->nlSys().copyOldSolutions();
247  _problem.getDisplacedProblem()->auxSys().copyOldSolutions();
248  }
249 
250  Real k_old = k;
252 
253  preIteration();
254  _problem.solve();
255  postIteration();
256 
257  // save the initial residual
258  if (iter == 0)
260 
261  // update eigenvalue
263  _eigenvalue = k;
264 
265  if (echo)
266  {
267  // output on screen the convergence history only when we want to and MOOSE output system is
268  // not used
269  keff_history.push_back(k);
270  if (solution_diff)
271  diff_history.push_back(*solution_diff);
272 
273  std::stringstream ss;
274  if (solution_diff)
275  {
276  ss << '\n';
277  ss << " +================+=====================+=====================+\n";
278  ss << " | iteration | eigenvalue | solution_difference |\n";
279  ss << " +================+=====================+=====================+\n";
280  unsigned int j = 0;
281  if (keff_history.size() > 10)
282  {
283  ss << " : : : :\n";
284  j = keff_history.size() - 10;
285  }
286  for (; j < keff_history.size(); j++)
287  ss << " | " << std::setw(14) << j << " | " << std::setw(19) << std::scientific
288  << std::setprecision(8) << keff_history[j] << " | " << std::setw(19) << std::scientific
289  << std::setprecision(8) << diff_history[j] << " |\n";
290  ss << " +================+=====================+=====================+\n" << std::flush;
291  }
292  else
293  {
294  ss << '\n';
295  ss << " +================+=====================+\n";
296  ss << " | iteration | eigenvalue |\n";
297  ss << " +================+=====================+\n";
298  unsigned int j = 0;
299  if (keff_history.size() > 10)
300  {
301  ss << " : : :\n";
302  j = keff_history.size() - 10;
303  }
304  for (; j < keff_history.size(); j++)
305  ss << " | " << std::setw(14) << j << " | " << std::setw(19) << std::scientific
306  << std::setprecision(8) << keff_history[j] << " |\n";
307  ss << " +================+=====================+\n" << std::flush;
308  ss << std::endl;
309  }
310  _console << ss.str();
311  }
312 
313  // increment iteration number here
314  iter++;
315 
316  if (cheb_on)
317  {
318  chebyshev(chebyshev_parameters, iter, solution_diff);
319  if (echo)
320  _console << " Chebyshev step: " << chebyshev_parameters.icheb << std::endl;
321  }
322 
323  if (echo)
324  _console
325  << " ________________________________________________________________________________ "
326  << std::endl;
327 
328  // not perform any convergence check when number of iterations is less than min_iter
329  if (iter >= min_iter)
330  {
331  // no need to check convergence of the last iteration
332  if (iter != max_iter)
333  {
334  bool converged = true;
335  Real keff_error = fabs(k_old - k) / k;
336  if (keff_error > tol_eig)
337  converged = false;
338  if (solution_diff)
339  if (*solution_diff > tol_x)
340  converged = false;
341  if (converged)
342  break;
343  }
344  else
345  break;
346  }
347  }
348 
349  // restore parameters changed by the executioner
350  _problem.es().parameters.set<Real>("linear solver tolerance") = tol1;
351  _problem.es().parameters.set<unsigned int>("nonlinear solver maximum iterations") = num1;
352  _problem.es().parameters.set<Real>("nonlinear solver relative residual tolerance") = tol2;
353 
354  // FIXME: currently power iteration use old and older solutions, so restore them
356  if (_problem.getDisplacedProblem() != NULL)
357  _problem.getDisplacedProblem()->restoreOldSolutions();
358 }
359 
360 void
362 {
363 }
364 
365 void
367 {
368 }
369 
370 void
372 {
373  if (getParam<bool>("output_before_normalization"))
374  {
375  _problem.timeStep()++;
376  Real t = _problem.time();
379  _problem.time() = t;
380  }
381 
382  Real s = 1.0;
384  {
385  _console << " Cannot let the normalization postprocessor on custom.\n";
386  _console << " Normalization is abandoned!" << std::endl;
387  }
388  else
389  {
391  if (!MooseUtils::absoluteFuzzyEqual(s, 1.0))
392  _console << " Solution is rescaled with factor " << s << " for normalization!" << std::endl;
393  }
394 
395  if ((!getParam<bool>("output_before_normalization")) || !MooseUtils::absoluteFuzzyEqual(s, 1.0))
396  {
397  _problem.timeStep()++;
398  Real t = _problem.time();
401  _problem.time() = t;
402  }
403 }
404 
405 Real
407 {
408  if (force)
410 
411  Real factor;
412  if (isParamValid("normal_factor"))
413  factor = getParam<Real>("normal_factor");
414  else
415  factor = _eigenvalue;
416  Real scaling = factor / _normalization;
417 
418  if (!MooseUtils::absoluteFuzzyEqual(scaling, 1.0))
419  {
420  // FIXME: we assume linear scaling here!
422  // update all aux variables and user objects
423  for (unsigned int i = 0; i < Moose::exec_types.size(); i++)
424  {
425  // EXEC_CUSTOM is special, should be treated only by specifically designed executioners.
426  if (Moose::exec_types[i] == EXEC_CUSTOM)
427  continue;
428 
430  }
431  }
432  return scaling;
433 }
434 
435 void
437 {
438  std::ostringstream ss;
439  ss << '\n';
440  ss << "*******************************************************\n";
441  ss << " Eigenvalue = " << std::fixed << std::setprecision(10) << _eigenvalue << '\n';
442  ss << "*******************************************************";
443 
444  _console << ss.str() << std::endl;
445 }
446 
448  : n_iter(50), fsmooth(2), finit(6), lgac(0), icheb(0), flux_error_norm_old(1), icho(0)
449 {
450 }
451 
452 void
454 {
455  finit = 6;
456  lgac = 0;
457  icheb = 0;
459  icho = 0;
460 }
461 
462 void
464  unsigned int iter,
465  const PostprocessorValue * solution_diff)
466 {
467  if (!solution_diff)
468  mooseError("solution diff is required for Chebyshev acceleration");
469 
470  if (chebyshev_parameters.lgac == 0)
471  {
472  if (chebyshev_parameters.icho == 0)
473  chebyshev_parameters.ratio = *solution_diff / chebyshev_parameters.flux_error_norm_old;
474  else
475  {
476  chebyshev_parameters.ratio = chebyshev_parameters.ratio_new;
477  chebyshev_parameters.icho = 0;
478  }
479 
480  if (iter > chebyshev_parameters.finit && chebyshev_parameters.ratio >= 0.4 &&
481  chebyshev_parameters.ratio <= 1)
482  {
483  chebyshev_parameters.lgac = 1;
484  chebyshev_parameters.icheb = 1;
485  chebyshev_parameters.error_begin = *solution_diff;
486  chebyshev_parameters.iter_begin = iter;
487  double alp = 2 / (2 - chebyshev_parameters.ratio);
488  std::vector<double> coef(2);
489  coef[0] = alp;
490  coef[1] = 1 - alp;
494  }
495  }
496  else
497  {
498  chebyshev_parameters.icheb++;
499  double gamma = acosh(2 / chebyshev_parameters.ratio - 1);
500  double alp = 4 / chebyshev_parameters.ratio *
501  std::cosh((chebyshev_parameters.icheb - 1) * gamma) /
502  std::cosh(chebyshev_parameters.icheb * gamma);
503  double beta = (1 - chebyshev_parameters.ratio / 2) * alp - 1;
504  /* if (iter<int(chebyshev_parameters.iter_begin+chebyshev_parameters.n_iter))
505  {
506  std::vector<double> coef(3);
507  coef[0] = alp;
508  coef[1] = 1-alp+beta;
509  coef[2] = -beta;
510  _eigen_sys.combineSystemSolution(NonlinearSystem::EIGEN, coef);
511  }
512  else
513  {*/
514  double gamma_new =
515  (*solution_diff / chebyshev_parameters.error_begin) *
516  (std::cosh((chebyshev_parameters.icheb - 1) * acosh(2 / chebyshev_parameters.ratio - 1)));
517  if (gamma_new < 1.0)
518  gamma_new = 1.0;
519 
520  chebyshev_parameters.ratio_new =
521  chebyshev_parameters.ratio / 2 *
522  (std::cosh(acosh(gamma_new) / (chebyshev_parameters.icheb - 1)) + 1);
523  if (gamma_new > 1.01)
524  {
525  chebyshev_parameters.lgac = 0;
526  // chebyshev_parameters.icheb = 0;
527  // if (chebyshev_parameters.icheb>30)
528  // {
529  if (chebyshev_parameters.icheb > 0)
530  {
531  chebyshev_parameters.icho = 1;
532  chebyshev_parameters.finit = iter;
533  }
534  else
535  {
536  chebyshev_parameters.icho = 0;
537  chebyshev_parameters.finit = iter + chebyshev_parameters.fsmooth;
538  }
539  }
540  else
541  {
542  std::vector<double> coef(3);
543  coef[0] = alp;
544  coef[1] = 1 - alp + beta;
545  coef[2] = -beta;
549  }
550  // }
551  }
552  chebyshev_parameters.flux_error_norm_old = *solution_diff;
553 }
554 
555 void
556 EigenExecutionerBase::nonlinearSolve(Real rel_tol, Real abs_tol, Real pfactor, Real & k)
557 {
558  makeBXConsistent(k);
559 
560  // turn on nonlinear flag so that eigen kernels opterate on the current solutions
562 
563  // set nonlinear solver controls
564  Real tol1 = _problem.es().parameters.get<Real>("nonlinear solver absolute residual tolerance");
565  Real tol2 = _problem.es().parameters.get<Real>("linear solver tolerance");
566  Real tol3 = _problem.es().parameters.get<Real>("nonlinear solver relative residual tolerance");
567 
568  _problem.es().parameters.set<Real>("nonlinear solver absolute residual tolerance") = abs_tol;
569  _problem.es().parameters.set<Real>("nonlinear solver relative residual tolerance") = rel_tol;
570  _problem.es().parameters.set<Real>("linear solver tolerance") = pfactor;
571 
572  // call nonlinear solve
573  _problem.solve();
574 
575  k = _source_integral;
576  _eigenvalue = k;
577 
578  _problem.es().parameters.set<Real>("nonlinear solver absolute residual tolerance") = tol1;
579  _problem.es().parameters.set<Real>("linear solver tolerance") = tol2;
580  _problem.es().parameters.set<Real>("nonlinear solver relative residual tolerance") = tol3;
581 }
Object is evaluated in every residual computation.
Definition: MooseTypes.h:96
void scaleSystemSolution(SYSTEMTAG tag, Real scaling_factor)
Scale the solution vector.
virtual Real & timeOld() const
EigenExecutionerBase(const InputParameters &parameters)
Constructor.
virtual void copyOldSolutions()
Shifts the solutions backwards in time.
Definition: SystemBase.C:693
NonlinearSystemBase & getNonlinearSystemBase()
For use with custom executioners that want to fire objects at a specific time.
Definition: MooseTypes.h:110
void addPrivateParam(const std::string &name, const T &value)
These method add a parameter to the InputParameters object which can be retrieved like any other para...
Object is evaluated only once at the beginning of the simulation.
Definition: MooseTypes.h:94
virtual void makeBXConsistent(Real k)
Normalize solution so that |Bx| = k.
const PostprocessorValue & getPostprocessorValueByName(const PostprocessorName &name)
Retrieve the value of the Postprocessor.
const Real & eigenvalueOld()
The old eigenvalue used by inverse power iterations.
virtual Real & time() const
The main MOOSE class responsible for handling user-defined parameters in almost every MOOSE system...
virtual void onTimestepEnd() override
virtual void postIteration()
Override this for actions that should take place after linear solve of each inverse power iteration...
virtual void checkIntegrity()
Make sure time kernel is not presented.
Object is evaluated at the end of every time step.
Definition: MooseTypes.h:100
void addRequiredParam(const std::string &name, const std::string &doc_string)
This method adds a parameter and documentation string to the InputParameters object that will be extr...
const PostprocessorValue & getPostprocessorValue(const std::string &name)
Retrieve the value of a Postprocessor or one of it&#39;s old or older values.
virtual EquationSystems & es() override
void buildSystemDoFIndices(SYSTEMTAG tag=ALL)
Build DoF indices for a system.
ExecFlagType execBitFlags() const
Build and return the execution flags as a bitfield.
void initSystemSolution(SYSTEMTAG tag, Real v)
Initialize the solution vector with a constant value.
void execute(const ExecFlagType &exec_type)
Convenience function for performing execution of MOOSE systems.
void eigenKernelOnOld()
Ask eigenkernels to operate on old or current solution vectors.
void combineSystemSolution(SYSTEMTAG tag, const std::vector< Real > &coefficients)
Linear combination of the solution vectors.
virtual std::shared_ptr< DisplacedProblem > getDisplacedProblem()
void setStartTime(const Real time)
Set the starting time for the simulation.
Definition: MooseApp.C:827
virtual void inversePowerIteration(unsigned int min_iter, unsigned int max_iter, Real pfactor, bool cheb_on, Real tol_eig, bool echo, PostprocessorName xdiff, Real tol_x, Real &k, Real &initial_res)
Perform inverse power iterations with the initial guess of the solution.
virtual void nonlinearSolve(Real rel_tol, Real abs_tol, Real pfactor, Real &k)
Perform nonlinear solve with the initial guess of the solution.
Real PostprocessorValue
MOOSE typedefs.
Definition: MooseTypes.h:73
virtual int & timeStep() const
InputParameters validParams< Executioner >()
Definition: Executioner.C:30
Executioners are objects that do the actual work of solving your problem.
Definition: Executioner.h:41
virtual Real & dt() const
bool isParamValid(const std::string &name) const
Test if the supplied parameter is valid.
Definition: MooseObject.h:67
AuxiliarySystem & getAuxiliarySystem()
const Real & _normalization
Postprocessor for normalization.
void chebyshev(Chebyshev_Parameters &params, unsigned int iter, const PostprocessorValue *solution_diff)
Real & _eigenvalue
Storage for the eigenvalue computed by the executioner.
bool containsEigenKernel() const
Weather or not the system contains eigen kernels.
virtual void restoreOldSolutions()
Restore old solutions from the backup vectors and deallocate them.
MooseEigenSystem & _eigen_sys
MooseApp & _app
The MooseApp this object is associated with.
Definition: MooseObject.h:108
InputParameters validParams< EigenExecutionerBase >()
virtual void printEigenvalue()
Print eigenvalue.
const InputParameters & parameters() const
Get the parameters of the object.
Definition: MooseObject.h:53
virtual void transient(bool trans)
virtual void saveOldSolutions()
Allocate vectors and save old solutions into them.
bool absoluteFuzzyEqual(const libMesh::Real &var1, const libMesh::Real &var2, const libMesh::Real &tol=libMesh::TOLERANCE *libMesh::TOLERANCE)
Function to check whether two variables are equal within an absolute tolerance.
void addParam(const std::string &name, const S &value, const std::string &doc_string)
These methods add an option parameter and a documentation string to the InputParameters object...
void mooseError(Args &&...args) const
Definition: MooseObject.h:80
const std::vector< ExecFlagType > exec_types
A static list of all the exec types.
Definition: MooseTypes.C:40
virtual Real normalizeSolution(bool force=true)
Normalize the solution vector based on the postprocessor value for normalization. ...
virtual void initialSetup()
ExecFlagType
Execution flags - when is the object executed/evaluated.
Definition: MooseTypes.h:90
const T & getUserObject(const std::string &name, unsigned int tid=0)
Get the user object by its name.
const ConsoleStream _console
An instance of helper class to write streams to the Console objects.
void initSystemSolutionOld(SYSTEMTAG tag, Real v)
virtual void solve() override
virtual void postExecute() override
Override this for actions that should take place after the main solve.
virtual void init() override
Initialize the executioner.
Base class for user-specific data.
Definition: UserObject.h:42
virtual void preIteration()
Override this for actions that should take place before linear solve of each inverse power iteration...
virtual void outputStep(ExecFlagType type)
Output the current step.
FEProblemBase & _fe_problem
Definition: Executioner.h:126
void addParamNamesToGroup(const std::string &space_delim_names, const std::string group_name)
This method takes a space delimited list of parameter names and adds them to the specified group name...
T & declareRestartableData(std::string data_name)
Declare a piece of data as "restartable".
Definition: Restartable.h:224