www.mooseframework.org
LStableDirk3.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 "LStableDirk3.h"
16 #include "NonlinearSystem.h"
17 #include "FEProblem.h"
18 #include "PetscSupport.h"
19 
20 template <>
23 {
25  return params;
26 }
27 
29  : TimeIntegrator(parameters),
30  _stage(1),
31  _gamma(-std::sqrt(2.) * std::cos(std::atan(std::sqrt(2.) / 4.) / 3.) / 2. +
32  std::sqrt(6.) * std::sin(std::atan(std::sqrt(2.) / 4.) / 3.) / 2. + 1.)
33 {
34  // Name the stage residuals "residual_stage1", "residual_stage2", etc.
35  for (unsigned int stage = 0; stage < 3; ++stage)
36  {
37  std::ostringstream oss;
38  oss << "residual_stage" << stage + 1;
39  _stage_residuals[stage] = &(_nl.addVector(oss.str(), false, GHOSTED));
40  }
41 
42  // Initialize parameters
43  _c[0] = _gamma;
44  _c[1] = .5 * (1 + _gamma);
45  _c[2] = 1.0;
46 
47  _a[0][0] = _gamma;
48  _a[1][0] = .5 * (1 - _gamma);
49  _a[1][1] = _gamma;
50  _a[2][0] = .25 * (-6 * _gamma * _gamma + 16 * _gamma - 1);
51  _a[2][1] = .25 * (6 * _gamma * _gamma - 20 * _gamma + 5);
52  _a[2][2] = _gamma;
53 }
54 
56 
57 void
59 {
60  // We are multiplying by the method coefficients in postStep(), so
61  // the time derivatives are of the same form at every stage although
62  // the current solution varies depending on the stage.
63  _u_dot = *_solution;
65  _u_dot *= 1. / _dt;
66  _u_dot.close();
67  _du_dot_du = 1. / _dt;
68 }
69 
70 void
72 {
73  // Time at end of step
74  Real time_old = _fe_problem.timeOld();
75 
76  // A for-loop would increment _stage too far, so we use an extra
77  // loop counter.
78  for (unsigned int current_stage = 1; current_stage < 4; ++current_stage)
79  {
80  // Set the current stage value
81  _stage = current_stage;
82 
83  // This ensures that all the Output objects in the OutputWarehouse
84  // have had solveSetup() called, and sets the default solver
85  // parameters for PETSc.
87 
88  _console << "Stage " << _stage << "\n";
89 
90  // Set the time for this stage
91  _fe_problem.time() = time_old + _c[_stage - 1] * _dt;
92 
93  // Do the solve
95  }
96 }
97 
98 void
99 LStableDirk3::postStep(NumericVector<Number> & residual)
100 {
101  // Error if _stage got messed up somehow.
102  if (_stage > 3)
103  mooseError("LStableDirk3::postStep(): Member variable _stage can only have values 1, 2, or 3.");
104 
105  // In the standard RK notation, the residual of stage 1 of s is given by:
106  //
107  // R := M*(Y_i - y_n)/dt - \sum_{j=1}^s a_{ij} * f(t_n + c_j*dt, Y_j) = 0
108  //
109  // where:
110  // .) M is the mass matrix
111  // .) Y_i is the stage solution
112  // .) dt is the timestep, and is accounted for in the _Re_time residual.
113  // .) f are the "non-time" residuals evaluated for a given stage solution.
114  // .) The minus signs are already "baked in" to the residuals and so do not appear below.
115 
116  // Store this stage's non-time residual. We are calling operator=
117  // here, and that calls close().
119 
120  // Build up the residual for this stage.
121  residual.add(1., _Re_time);
122  for (unsigned int j = 0; j < _stage; ++j)
123  residual.add(_a[_stage - 1][j], *_stage_residuals[j]);
124  residual.close();
125 }
NonlinearSystemBase & _nl
virtual Real & timeOld() const
void initPetscOutput()
Reinitialize petsc output for proper linear/nonlinear iteration display.
virtual void solve()
Definition: LStableDirk3.C:71
NumericVector< Number > & _u_dot
solution vector for u^dot
NonlinearSystemBase & getNonlinearSystemBase()
virtual NumericVector< Number > & addVector(const std::string &vector_name, const bool project, const ParallelType type)
Adds a solution length vector to the system.
Definition: SystemBase.C:490
FEProblemBase & _fe_problem
NumericVector< Number > * _stage_residuals[3]
Definition: LStableDirk3.h:67
NumericVector< Number > & _Re_non_time
residual vector for non-time contributions
virtual Real & time() const
The main MOOSE class responsible for handling user-defined parameters in almost every MOOSE system...
InputParameters validParams< LStableDirk3 >()
Definition: LStableDirk3.C:22
Real _a[3][3]
Definition: LStableDirk3.h:81
InputParameters validParams< TimeIntegrator >()
virtual void computeTimeDerivatives()
Definition: LStableDirk3.C:58
unsigned int _stage
Definition: LStableDirk3.h:64
virtual ~LStableDirk3()
Definition: LStableDirk3.C:55
Real & _du_dot_du
solution vector for
Base class for time integrators.
virtual void postStep(NumericVector< Number > &residual)
Callback to the TimeIntegrator called immediately after the residuals are computed in NonlinearSystem...
Definition: LStableDirk3.C:99
const Real _gamma
Definition: LStableDirk3.h:70
virtual System & system() override
Get the reference to the libMesh system.
LStableDirk3(const InputParameters &parameters)
Definition: LStableDirk3.C:28
NumericVector< Number > & _Re_time
residual vector for time contributions
void mooseError(Args &&...args) const
Definition: MooseObject.h:80
const NumericVector< Number > & _solution_old
Real _c[3]
Definition: LStableDirk3.h:74
const ConsoleStream _console
An instance of helper class to write streams to the Console objects.
const NumericVector< Number > *& _solution
solution vectors