Skip to content

Commit

Permalink
Add explicit Heun, Ralston, and Midpoint time integration methods.
Browse files Browse the repository at this point in the history 
* These three methods have a common ExplicitRK2 base class.
* The original RungeKutta2 method is equivalent to ExplicitMidpoint.
* Improve documentation.

Refs idaholab#1929.
  • Loading branch information
@jwpeterson
jwpeterson committed Sep 15, 2015
1 parent bc5cfc1 commit 78842f5
Show file tree
Hide file tree
Showing 10 changed files with 508 additions and 0 deletions.
50 changes: 50 additions & 0 deletions framework/include/timeintegrators/ExplicitMidpoint.h
View file
@@ -0,0 +1,50 @@
/****************************************************************/
/* DO NOT MODIFY THIS HEADER */
/* MOOSE - Multiphysics Object Oriented Simulation Environment */
/* */
/* (c) 2010 Battelle Energy Alliance, LLC */
/* ALL RIGHTS RESERVED */
/* */
/* Prepared by Battelle Energy Alliance, LLC */
/* Under Contract No. DE-AC07-05ID14517 */
/* With the U. S. Department of Energy */
/* */
/* See COPYRIGHT for full restrictions */
/****************************************************************/

#ifndef EXPLICITMIDPOINT_H
#define EXPLICITMIDPOINT_H

#include "ExplicitRK2.h"

class ExplicitMidpoint;

template<>
InputParameters validParams<ExplicitMidpoint>();

/**
* The explicit midpoint time integration method.
*
* The Butcher tableau for this method is:
* 0 | 0
* 1/2 | 1/2 0
* ---------------------
* | 0 1
*
* See: ExplicitRK2.h for more information.
*/
class ExplicitMidpoint : public ExplicitRK2
{
public:
ExplicitMidpoint(const InputParameters & parameters);
virtual ~ExplicitMidpoint() {}

protected:
/// Method coefficient overrides
virtual Real a() const { return .5; }
virtual Real b1() const { return 0.; }
virtual Real b2() const { return 1.; }
};


#endif /* EXPLICITMIDPOINT_H */
92 changes: 92 additions & 0 deletions framework/include/timeintegrators/ExplicitRK2.h
View file
@@ -0,0 +1,92 @@
/****************************************************************/
/* DO NOT MODIFY THIS HEADER */
/* MOOSE - Multiphysics Object Oriented Simulation Environment */
/* */
/* (c) 2010 Battelle Energy Alliance, LLC */
/* ALL RIGHTS RESERVED */
/* */
/* Prepared by Battelle Energy Alliance, LLC */
/* Under Contract No. DE-AC07-05ID14517 */
/* With the U. S. Department of Energy */
/* */
/* See COPYRIGHT for full restrictions */
/****************************************************************/

#ifndef EXPLICITRK2_H
#define EXPLICITRK2_H

#include "TimeIntegrator.h"

class ExplicitRK2;

template<>
InputParameters validParams<ExplicitRK2>();

/**
* Base class for three different explicit second-order Runge-Kutta
* time integration methods:
* - Explicit midpoint method (ExplicitMidpoint.C)
* - Heun's method (Heun.C)
* - Ralston's method (Ralston.C)
*
* Each of these methods is characterized by the following generic
* Butcher tableau:
* 0 | 0
* a | a 0
* ---------------------
* | b1 b2
*
* where a, b1, and b2 are constants that define the different
* methods.
*
* The stability function for each of these methods is:
* R(z) = z*(z + 2)/2 + 1
*
* The methods are all explicit, so they are neither A-stable nor L-stable,
* lim R(z), z->oo = oo
*
* See also:
* https://en.wikipedia.org/wiki/Runge%E2%80%93Kutta_methods
* https://en.wikipedia.org/wiki/Midpoint_method
* https://en.wikipedia.org/wiki/Heun%27s_method
*
* All three methods require two mass matrix (linear) system solves
* per timestep. Although strictly speaking these are "two stage"
* methods, we treat the "update" step as a third stage, since in
* finite element analysis the update step requires a mass matrix
* solve.
*
* IMPORTANT: To use the explicit TimeIntegrators derived from this
* method, you must generally add "implicit=false" to the Kernels,
* Materials, etc. used in your simulation, so that MOOSE evaluates
* them correctly! An important exception are TimeDerivative kernels,
* which should never be marked "implicit=false".
*/
class ExplicitRK2 : public TimeIntegrator
{
public:
ExplicitRK2(const InputParameters & parameters);
virtual ~ExplicitRK2();

virtual int order() { return 2; }

virtual void computeTimeDerivatives();
virtual void solve();
virtual void postStep(NumericVector<Number> & residual);

protected:
unsigned int _stage;

/// Buffer to store non-time residual from the first stage.
NumericVector<Number> & _residual_old;

/// The method coefficients. See the table above for description.
/// These are pure virtual in the base class, and must be overridden
/// in subclasses to implement different schemes.
virtual Real a() const = 0;
virtual Real b1() const = 0;
virtual Real b2() const = 0;
};


#endif /* EXPLICITRK2_H */
50 changes: 50 additions & 0 deletions framework/include/timeintegrators/Heun.h
View file
@@ -0,0 +1,50 @@
/****************************************************************/
/* DO NOT MODIFY THIS HEADER */
/* MOOSE - Multiphysics Object Oriented Simulation Environment */
/* */
/* (c) 2010 Battelle Energy Alliance, LLC */
/* ALL RIGHTS RESERVED */
/* */
/* Prepared by Battelle Energy Alliance, LLC */
/* Under Contract No. DE-AC07-05ID14517 */
/* With the U. S. Department of Energy */
/* */
/* See COPYRIGHT for full restrictions */
/****************************************************************/

#ifndef HEUN_H
#define HEUN_H

#include "ExplicitRK2.h"

class Heun;

template<>
InputParameters validParams<Heun>();

/**
* Heun's (aka improved Euler) time integration method.
*
* The Butcher tableau for this method is:
* 0 | 0
* 1 | 1 0
* ---------------------
* | 1/2 1/2
*
* See: ExplicitRK2.h for more information.
*/
class Heun : public ExplicitRK2
{
public:
Heun(const InputParameters & parameters);
virtual ~Heun() {}

protected:
/// Method coefficient overrides
virtual Real a() const { return 1.; }
virtual Real b1() const { return .5; }
virtual Real b2() const { return .5; }
};


#endif /* HEUN_H */
50 changes: 50 additions & 0 deletions framework/include/timeintegrators/Ralston.h
View file
@@ -0,0 +1,50 @@
/****************************************************************/
/* DO NOT MODIFY THIS HEADER */
/* MOOSE - Multiphysics Object Oriented Simulation Environment */
/* */
/* (c) 2010 Battelle Energy Alliance, LLC */
/* ALL RIGHTS RESERVED */
/* */
/* Prepared by Battelle Energy Alliance, LLC */
/* Under Contract No. DE-AC07-05ID14517 */
/* With the U. S. Department of Energy */
/* */
/* See COPYRIGHT for full restrictions */
/****************************************************************/

#ifndef RALSTON_H
#define RALSTON_H

#include "ExplicitRK2.h"

class Ralston;

template<>
InputParameters validParams<Ralston>();

/**
* Ralston's time integration method.
*
* The Butcher tableau for this method is:
* 0 | 0
* 2/3 | 2/3 0
* ---------------------
* | 1/4 3/4
*
* See: ExplicitRK2.h for more information.
*/
class Ralston : public ExplicitRK2
{
public:
Ralston(const InputParameters & parameters);
virtual ~Ralston() {}

protected:
/// Method coefficient overrides
virtual Real a() const { return 2./3.; }
virtual Real b1() const { return .25; }
virtual Real b2() const { return .75; }
};


#endif /* RALSTON_H */
17 changes: 17 additions & 0 deletions framework/include/timeintegrators/TimeIntegrator.h
View file
Expand Up @@ -48,7 +48,24 @@ class TimeIntegrator :
virtual void preSolve() { }
virtual void preStep() { }
virtual void solve();

/**
* Callback to the TimeIntegrator called immediately after the
* residuals are computed in NonlinearSystem::computeResidual() (it
* is not really named well...). The residual vector which is
* passed in to this function should be filled in by the user with
* the _Re_time and _Re_non_time vectors in a way that makes sense
* for the particular TimeIntegration method.
*/
virtual void postStep(NumericVector<Number> & /*residual*/) { }

/**
* Callback to the TimeIntegrator called immediately after
* TimeIntegrator::solve() (so the name does make sense!). See
* e.g. CrankNicolson for an example of what can be done in the
* postSolve() callback -- there it is used to move the residual
* from the "old" timestep forward in time to avoid recomputing it.
*/
virtual void postSolve() { }

virtual int order() = 0;
Expand Down
6 changes: 6 additions & 0 deletions framework/src/base/Moose.C
View file
Expand Up @@ -287,9 +287,12 @@
#include "CrankNicolson.h"
#include "ExplicitEuler.h"
#include "RungeKutta2.h"
#include "ExplicitMidpoint.h"
#include "Dirk.h"
#include "LStableDirk2.h"
#include "ImplicitMidpoint.h"
#include "Heun.h"
#include "Ralston.h"
//
#include "SimplePredictor.h"
#include "AdamsPredictor.h"
Expand Down Expand Up @@ -675,9 +678,12 @@ registerObjects(Factory & factory)
registerTimeIntegrator(CrankNicolson);
registerTimeIntegrator(ExplicitEuler);
registerTimeIntegrator(RungeKutta2);
registerTimeIntegrator(ExplicitMidpoint);
registerDeprecatedObjectName(Dirk, "Dirk", "09/22/2015 12:00");
registerTimeIntegrator(LStableDirk2);
registerTimeIntegrator(ImplicitMidpoint);
registerTimeIntegrator(Heun);
registerTimeIntegrator(Ralston);
// predictors
registerPredictor(SimplePredictor);
registerPredictor(AdamsPredictor);
Expand Down
28 changes: 28 additions & 0 deletions framework/src/timeintegrators/ExplicitMidpoint.C
View file
@@ -0,0 +1,28 @@
/****************************************************************/
/* DO NOT MODIFY THIS HEADER */
/* MOOSE - Multiphysics Object Oriented Simulation Environment */
/* */
/* (c) 2010 Battelle Energy Alliance, LLC */
/* ALL RIGHTS RESERVED */
/* */
/* Prepared by Battelle Energy Alliance, LLC */
/* Under Contract No. DE-AC07-05ID14517 */
/* With the U. S. Department of Energy */
/* */
/* See COPYRIGHT for full restrictions */
/****************************************************************/

#include "ExplicitMidpoint.h"

template<>
InputParameters validParams<ExplicitMidpoint>()
{
InputParameters params = validParams<ExplicitRK2>();

return params;
}

ExplicitMidpoint::ExplicitMidpoint(const InputParameters & parameters) :
ExplicitRK2(parameters)
{
}

0 comments on commit 78842f5

Please sign in to comment.