Commit 5957f906 authored by Dmitry Yu. Naumov's avatar Dmitry Yu. Naumov
Browse files

[NL] Remove unused time discretization schemes.

This are all but backward Euler.
The few testcases in H2-pp were updated to backward Euler
without loss of accuracy.

 - All processes using a Newton scheme are discretized
for the backward Euler scheme only. They are not usable
with other time discretizations.
 - Greatly simplifies code maintenance in this part.
parent 9ce7bf3e
......@@ -44,13 +44,11 @@ void MatrixTranslatorGeneral<ODESystemTag::FirstOrderImplicitQuasilinear>::
void MatrixTranslatorGeneral<ODESystemTag::FirstOrderImplicitQuasilinear>::
computeResidual(GlobalMatrix const& M, GlobalMatrix const& K,
GlobalVector const& b, GlobalVector const& x_new_timestep,
GlobalVector const& b, GlobalVector const& x_curr,
GlobalVector const& xdot, GlobalVector& res) const
{
namespace LinAlg = MathLib::LinAlg;
auto const& x_curr = _time_disc.getCurrentX(x_new_timestep);
// res = M * x_dot + K * x_curr - b
LinAlg::matMult(M, xdot, res); // the local vector x_dot seems to be
// necessary because of this multiplication
......@@ -66,161 +64,4 @@ void MatrixTranslatorGeneral<ODESystemTag::FirstOrderImplicitQuasilinear>::
LinAlg::copy(Jac_in, Jac_out);
}
void MatrixTranslatorForwardEuler<ODESystemTag::FirstOrderImplicitQuasilinear>::
computeA(GlobalMatrix const& M, GlobalMatrix const& /*K*/,
GlobalMatrix& A) const
{
namespace LinAlg = MathLib::LinAlg;
auto const dxdot_dx = _fwd_euler.getNewXWeight();
// A = M * dxdot_dx
LinAlg::copy(M, A);
LinAlg::scale(A, dxdot_dx);
}
void MatrixTranslatorForwardEuler<ODESystemTag::FirstOrderImplicitQuasilinear>::
computeRhs(const GlobalMatrix& M, const GlobalMatrix& K,
const GlobalVector& b, GlobalVector& rhs) const
{
namespace LinAlg = MathLib::LinAlg;
auto& tmp = NumLib::GlobalVectorProvider::provider.getVector(_tmp_id);
_fwd_euler.getWeightedOldX(tmp);
auto const& x_old = _fwd_euler.getXOld();
// rhs = b + M * weighted_old_x - K * x_old
LinAlg::matMult(K, x_old, rhs); // rhs = K * x_old
LinAlg::aypx(rhs, -1.0, b); // rhs = b - K * x_old
LinAlg::matMultAdd(M, tmp, rhs, rhs); // rhs += M * weighted_old_x
NumLib::GlobalVectorProvider::provider.releaseVector(tmp);
}
void MatrixTranslatorForwardEuler<ODESystemTag::FirstOrderImplicitQuasilinear>::
computeResidual(GlobalMatrix const& M, GlobalMatrix const& K,
GlobalVector const& b, GlobalVector const& x_new_timestep,
GlobalVector const& xdot, GlobalVector& res) const
{
namespace LinAlg = MathLib::LinAlg;
auto const& x_curr = _fwd_euler.getCurrentX(x_new_timestep);
// res = M * x_dot + K * x_curr - b
LinAlg::matMult(M, xdot, res);
LinAlg::matMultAdd(K, x_curr, res, res);
LinAlg::axpy(res, -1.0, b);
}
void MatrixTranslatorForwardEuler<ODESystemTag::FirstOrderImplicitQuasilinear>::
computeJacobian(GlobalMatrix const& Jac_in, GlobalMatrix& Jac_out) const
{
namespace LinAlg = MathLib::LinAlg;
LinAlg::copy(Jac_in, Jac_out);
}
void MatrixTranslatorCrankNicolson<
ODESystemTag::FirstOrderImplicitQuasilinear>::
computeA(GlobalMatrix const& M, GlobalMatrix const& K,
GlobalMatrix& A) const
{
namespace LinAlg = MathLib::LinAlg;
auto const dxdot_dx = _crank_nicolson.getNewXWeight();
auto const theta = _crank_nicolson.getTheta();
// A = theta * (M * dxdot_dx + K) + dxdot_dx * _M_bar
LinAlg::copy(M, A);
LinAlg::aypx(A, dxdot_dx, K); // A = M * dxdot_dx + K
LinAlg::scale(A, theta); // A *= theta
LinAlg::axpy(A, dxdot_dx, _M_bar); // A += dxdot_dx * _M_bar
}
void MatrixTranslatorCrankNicolson<
ODESystemTag::FirstOrderImplicitQuasilinear>::
computeRhs(const GlobalMatrix& M, const GlobalMatrix& /*K*/,
const GlobalVector& b, GlobalVector& rhs) const
{
namespace LinAlg = MathLib::LinAlg;
auto& tmp = NumLib::GlobalVectorProvider::provider.getVector(_tmp_id);
_crank_nicolson.getWeightedOldX(tmp);
auto const theta = _crank_nicolson.getTheta();
// rhs = theta * (b + M * weighted_old_x) + _M_bar * weighted_old_x -
// _b_bar;
LinAlg::matMultAdd(M, tmp, b, rhs); // rhs = b + M * weighted_old_x
LinAlg::scale(rhs, theta); // rhs *= theta
LinAlg::matMultAdd(_M_bar, tmp, rhs, rhs); // rhs += _M_bar * weighted_old_x
LinAlg::axpy(rhs, -1.0, _b_bar); // rhs -= b
NumLib::GlobalVectorProvider::provider.releaseVector(tmp);
}
void MatrixTranslatorCrankNicolson<
ODESystemTag::FirstOrderImplicitQuasilinear>::
computeResidual(GlobalMatrix const& M, GlobalMatrix const& K,
GlobalVector const& b, GlobalVector const& x_new_timestep,
GlobalVector const& xdot, GlobalVector& res) const
{
namespace LinAlg = MathLib::LinAlg;
auto const& x_curr = _crank_nicolson.getCurrentX(x_new_timestep);
auto const theta = _crank_nicolson.getTheta();
// res = theta * (M * x_dot + K*x_curr - b) + _M_bar * x_dot + _b_bar
LinAlg::matMult(M, xdot, res); // res = M * x_dot
LinAlg::matMultAdd(K, x_curr, res, res); // res += K * x_curr
LinAlg::axpy(res, -1.0, b); // res = M * x_dot + K * x_curr - b
LinAlg::aypx(res, theta, _b_bar); // res = res * theta + _b_bar
LinAlg::matMultAdd(_M_bar, xdot, res, res); // rs += _M_bar * x_dot
}
void MatrixTranslatorCrankNicolson<
ODESystemTag::FirstOrderImplicitQuasilinear>::
computeJacobian(GlobalMatrix const& Jac_in, GlobalMatrix& Jac_out) const
{
namespace LinAlg = MathLib::LinAlg;
auto const dxdot_dx = _crank_nicolson.getNewXWeight();
auto const theta = _crank_nicolson.getTheta();
// J = theta * Jac + dxdot_dx * _M_bar
LinAlg::copy(Jac_in, Jac_out);
LinAlg::scale(Jac_out, theta);
LinAlg::axpy(Jac_out, dxdot_dx, _M_bar);
}
void MatrixTranslatorCrankNicolson<
ODESystemTag::FirstOrderImplicitQuasilinear>::
pushMatrices(GlobalMatrix const& M, GlobalMatrix const& K,
GlobalVector const& b)
{
namespace LinAlg = MathLib::LinAlg;
auto const theta = _crank_nicolson.getTheta();
// Note: using x_old here is correct, since this method is called from
// within
// CrankNicolson::pushState() __after__ x_old has been updated to the
// result
// from the timestep just finished.
auto const& x_old = _crank_nicolson.getXOld();
// _M_bar = (1.0-theta) * M;
LinAlg::copy(M, _M_bar);
LinAlg::scale(_M_bar, 1.0 - theta);
// _b_bar = (1.0-theta) * (K * x_old - b)
LinAlg::matMult(K, x_old, _b_bar);
LinAlg::axpy(_b_bar, -1.0, b);
LinAlg::scale(_b_bar, 1.0 - theta);
}
} // namespace NumLib
......@@ -62,19 +62,6 @@ public:
virtual void computeJacobian(GlobalMatrix const& Jac_in,
GlobalMatrix& Jac_out) const = 0;
/*! Allows to store the given matrices internally for later use.
*
* \remark
* This method has been provided in order to be able to implement the
* CrankNicolson
* scheme.
*/
virtual void pushMatrices(GlobalMatrix const& /*M*/,
GlobalMatrix const& /*K*/,
GlobalVector const& /*b*/)
{
}
virtual ~MatrixTranslator() = default;
};
......@@ -136,172 +123,12 @@ private:
mutable std::size_t _tmp_id = 0u;
};
/*! GlobalMatrix translator used with the ForwardEuler scheme.
*
* \tparam ODETag a tag indicating the type of equation.
*/
template <ODESystemTag ODETag>
class MatrixTranslatorForwardEuler;
/*! GlobalMatrix translator for first order implicit quasi-linear ODEs,
* used with the ForwardEuler scheme.
*
* \see ODESystemTag::FirstOrderImplicitQuasilinear
* \remark
* You might also want read the remarks on
* \ref concept_time_discretization "time discretization".
*/
template <>
class MatrixTranslatorForwardEuler<ODESystemTag::FirstOrderImplicitQuasilinear>
: public MatrixTranslator<ODESystemTag::FirstOrderImplicitQuasilinear>
{
public:
/*! Constructs a new instance.
*
* \param timeDisc the time discretization scheme to be used.
*/
explicit MatrixTranslatorForwardEuler(ForwardEuler const& timeDisc)
: _fwd_euler(timeDisc)
{
}
//! Computes \f$ A = M \cdot \alpha \f$.
void computeA(GlobalMatrix const& M, GlobalMatrix const& /*K*/,
GlobalMatrix& A) const override;
//! Computes \f$ \mathtt{rhs} = M \cdot x_O - K \cdot x_O + b \f$.
void computeRhs(const GlobalMatrix& M, const GlobalMatrix& K,
const GlobalVector& b, GlobalVector& rhs) const override;
//! Computes \f$ r = M \cdot \hat x + K \cdot x_C - b \f$.
void computeResidual(GlobalMatrix const& M, GlobalMatrix const& K,
GlobalVector const& b,
GlobalVector const& x_new_timestep,
GlobalVector const& xdot,
GlobalVector& res) const override;
//! Writes \c Jac_in to \c Jac_out.
//! \todo Do not copy.
void computeJacobian(GlobalMatrix const& Jac_in,
GlobalMatrix& Jac_out) const override;
private:
ForwardEuler const& _fwd_euler; //!< the time discretization used.
//! ID of the vector storing intermediate computations.
mutable std::size_t _tmp_id = 0u;
};
/*! GlobalMatrix translator used with the CrankNicolson scheme.
*
* \tparam ODETag a tag indicating the type of equation.
*/
template <ODESystemTag ODETag>
class MatrixTranslatorCrankNicolson;
/*! GlobalMatrix translator for first order implicit quasi-linear ODEs,
* used with the CrankNicolson scheme.
*
* \see ODESystemTag::FirstOrderImplicitQuasilinear
* \remark
* You might also want read the remarks on
* \ref concept_time_discretization "time discretization".
*/
template <>
class MatrixTranslatorCrankNicolson<ODESystemTag::FirstOrderImplicitQuasilinear>
: public MatrixTranslator<ODESystemTag::FirstOrderImplicitQuasilinear>
{
public:
/*! Constructs a new instance.
*
* \param timeDisc the time discretization scheme to be used.
*/
explicit MatrixTranslatorCrankNicolson(CrankNicolson const& timeDisc)
: _crank_nicolson(timeDisc),
_M_bar(NumLib::GlobalMatrixProvider::provider.getMatrix()),
_b_bar(NumLib::GlobalVectorProvider::provider.getVector())
{
}
~MatrixTranslatorCrankNicolson() override
{
NumLib::GlobalMatrixProvider::provider.releaseMatrix(
_M_bar);
NumLib::GlobalVectorProvider::provider.releaseVector(
_b_bar);
}
//! Computes \f$ A = \theta \cdot (M \cdot \alpha + K) + \bar M \cdot \alpha
//! \f$.
void computeA(GlobalMatrix const& M, GlobalMatrix const& K,
GlobalMatrix& A) const override;
//! Computes \f$ \mathtt{rhs} = \theta \cdot (M \cdot x_O + b) + \bar M
//! \cdot x_O - \bar b \f$.
void computeRhs(const GlobalMatrix& M, const GlobalMatrix& /*K*/,
const GlobalVector& b, GlobalVector& rhs) const override;
//! Computes \f$ r = \theta \cdot (M \cdot \hat x + K \cdot x_C - b) + \bar
//! M \cdot \hat x + \bar b \f$.
void computeResidual(GlobalMatrix const& M, GlobalMatrix const& K,
GlobalVector const& b,
GlobalVector const& x_new_timestep,
GlobalVector const& xdot,
GlobalVector& res) const override;
/*! Computes \f$ \mathtt{Jac\_out} = \theta \cdot \mathtt{Jac\_in} + \bar M
* \cdot \alpha \f$.
*
* Where \c Jac_in is the Jacobian as assembled by the ODE system, i.e. in
* the same fashion as for the BackwardEuler scheme.
*/
void computeJacobian(GlobalMatrix const& Jac_in,
GlobalMatrix& Jac_out) const override;
/*! Saves internal state for use in the successive timestep;
* computes \f$ \bar M \f$ and \f$ \bar b \f$.
*
* \f$ \bar M \f$ and \f$ \bar b \f$ are computed as follows:
* \f{align}{
* \bar M &= (1-\theta) \cdot M \\
* \bar b &= (1-\theta) \cdot ( K \cdot x_n - b )
* \f}
*
* Where \f$ x_n \f$ is the solution at the timestep just finished.
*/
void pushMatrices(GlobalMatrix const& M, GlobalMatrix const& K,
GlobalVector const& b) override;
private:
CrankNicolson const& _crank_nicolson;
GlobalMatrix&
_M_bar; //!< Used to adjust matrices and vectors assembled by the ODE.
//!< \see pushMatrices()
GlobalVector& _b_bar; //!< Used to adjust vectors assembled by the ODE.
//!< \see pushMatrices()
//! ID of the vector storing intermediate computations.
mutable std::size_t _tmp_id = 0u;
};
//! Creates a GlobalMatrix translator suitable to work together with the given
//! time discretization scheme.
template <ODESystemTag ODETag>
std::unique_ptr<MatrixTranslator<ODETag>> createMatrixTranslator(
TimeDiscretization const& timeDisc)
{
if (auto* fwd_euler = dynamic_cast<ForwardEuler const*>(&timeDisc))
{
return std::unique_ptr<MatrixTranslator<ODETag>>(
new MatrixTranslatorForwardEuler<ODETag>(*fwd_euler));
}
if (auto* crank = dynamic_cast<CrankNicolson const*>(&timeDisc))
{
return std::unique_ptr<MatrixTranslator<ODETag>>(
new MatrixTranslatorCrankNicolson<ODETag>(*crank));
}
return std::unique_ptr<MatrixTranslator<ODETag>>(
new MatrixTranslatorGeneral<ODETag>(timeDisc));
}
......
......@@ -22,12 +22,6 @@
namespace NumLib
{
void NonlinearSolver<NonlinearSolverTag::Picard>::assemble(
std::vector<GlobalVector*> const& x, int const process_id) const
{
_equation_system->assemble(x, process_id);
}
void NonlinearSolver<NonlinearSolverTag::Picard>::
calculateNonEquilibriumInitialResiduum(std::vector<GlobalVector*> const& x,
int const process_id)
......@@ -211,15 +205,6 @@ NonlinearSolverStatus NonlinearSolver<NonlinearSolverTag::Picard>::solve(
return {error_norms_met, iteration};
}
void NonlinearSolver<NonlinearSolverTag::Newton>::assemble(
std::vector<GlobalVector*> const& x, int const process_id) const
{
_equation_system->assemble(x, process_id);
// TODO if the equation system would be reset to nullptr after each
// assemble() or solve() call, the user would be forced to set the
// equation every time and could not forget it.
}
void NonlinearSolver<NonlinearSolverTag::Newton>::
calculateNonEquilibriumInitialResiduum(std::vector<GlobalVector*> const& x,
int const process_id)
......
......@@ -35,19 +35,6 @@ namespace NumLib
class NonlinearSolverBase
{
public:
/*! Only assemble the equation system.
*
* \note This method is needed to preload CrankNicolson time discretization
* scheme. It is not used for the general solver steps; in those only the
* solve() method is needed.
*
* \param x the state at which the equation system will be
* assembled.
* \param process_id the process' id which will be assembled.
*/
virtual void assemble(std::vector<GlobalVector*> const& x,
int const process_id) const = 0;
virtual void calculateNonEquilibriumInitialResiduum(
std::vector<GlobalVector*> const& x, int const process_id) = 0;
......@@ -112,9 +99,6 @@ public:
_convergence_criterion = &conv_crit;
}
void assemble(std::vector<GlobalVector*> const& x,
int const process_id) const override;
void calculateNonEquilibriumInitialResiduum(
std::vector<GlobalVector*> const& x, int const process_id) override;
......@@ -190,9 +174,6 @@ public:
_convergence_criterion = &conv_crit;
}
void assemble(std::vector<GlobalVector*> const& x,
int const process_id) const override;
void calculateNonEquilibriumInitialResiduum(
std::vector<GlobalVector*> const& x, int const process_id) override;
......
......@@ -59,87 +59,4 @@ double BackwardEuler::getRelativeChangeFromPreviousTimestep(
return computeRelativeChangeFromPreviousTimestep(x, _x_old, norm_type);
}
double ForwardEuler::getRelativeChangeFromPreviousTimestep(
GlobalVector const& x, MathLib::VecNormType norm_type)
{
return computeRelativeChangeFromPreviousTimestep(x, _x_old, norm_type);
}
double CrankNicolson::getRelativeChangeFromPreviousTimestep(
GlobalVector const& x, MathLib::VecNormType norm_type)
{
return computeRelativeChangeFromPreviousTimestep(x, _x_old, norm_type);
}
double BackwardDifferentiationFormula::getRelativeChangeFromPreviousTimestep(
GlobalVector const& x, MathLib::VecNormType norm_type)
{
return computeRelativeChangeFromPreviousTimestep(
x, *_xs_old[_offset], norm_type);
}
void BackwardDifferentiationFormula::pushState(
const double /*t*/,
GlobalVector const& x,
InternalMatrixStorage const& /*strg*/)
{
namespace LinAlg = MathLib::LinAlg;
// TODO use boost circular buffer?
// until _xs_old is filled, lower-order BDF formulas are used.
if (_xs_old.size() < _num_steps)
{
_xs_old.push_back(&NumLib::GlobalVectorProvider::provider.getVector(x));
}
else
{
LinAlg::copy(x, *_xs_old[_offset]);
_offset =
(_offset + 1) % _num_steps; // treat _xs_old as a circular buffer
}
}
namespace detail
{
//! Coefficients used in the backward differentiation formulas.
const double BDF_Coeffs[6][7] = {
// leftmost column: weight of the solution at the new timestep
// signs of columns > 1 are flipped compared to standard BDF tableaus
{1.0, 1.0},
{1.5, 2.0, -0.5},
{11.0 / 6.0, 3.0, -1.5, 1.0 / 3.0},
{25.0 / 12.0, 4.0, -3.0, 4.0 / 3.0, -0.25},
{137.0 / 60.0, 5.0, -5.0, 10.0 / 3.0, -1.25, 1.0 / 5.0},
{147.0 / 60.0, 6.0, -7.5, 20.0 / 3.0, -3.75, 6.0 / 5.0, -1.0 / 6.0}
// coefficient of (for BDF(6), the oldest state, x_n, is always rightmost)
// x_+6, x_+5, x_+4, x_+3, x_+2, x_+1, x_n
};
} // namespace detail
double BackwardDifferentiationFormula::getNewXWeight() const
{
auto const k = eff_num_steps();
return detail::BDF_Coeffs[k - 1][0] / _delta_t;
}
void BackwardDifferentiationFormula::getWeightedOldX(GlobalVector& y) const
{
namespace LinAlg = MathLib::LinAlg;
auto const k = eff_num_steps();
auto const* const BDFk = detail::BDF_Coeffs[k - 1];
// compute linear combination \sum_{i=0}^{k-1} BDFk_{k-i} \cdot x_{n+i}
LinAlg::copy(*_xs_old[_offset], y); // _xs_old[offset] = x_n
LinAlg::scale(y, BDFk[k]);
for (unsigned i = 1; i < k; ++i)
{
auto const off = (_offset + i) % k;
LinAlg::axpy(y, BDFk[k - i], *_xs_old[off]);
}
LinAlg::scale(y, 1.0 / _delta_t);
}
} // end of namespace NumLib
......@@ -21,28 +21,6 @@ namespace NumLib
//! \addtogroup ODESolver
//! @{
//! Interface that allows managing an additional internal state as required by
//! certain time discretization schemes.
class InternalMatrixStorage
{
public:
/*! Triggers a refresh of the internal matrix/vector storage.
*
* \remark
* This method is needed in particular to fully implement the
* interaction of the CrankNicolson scheme with other classes.
*
* \attention
* This method must be called (if it is called) from within
* TimeDiscretization::pushState() \b after the internal state of
* the TimeDiscretization has been set to the new solution.
* Otherwise the pushMatrices() method of MatrixTranslator's will break!
*/
virtual void pushMatrices() const = 0;
virtual ~InternalMatrixStorage() = default;
};
/*! Interface of time discretization schemes for first-order ODEs.
*
* The purpose of TimeDiscretization instances is to store the solution history of
......@@ -103,13 +81,8 @@ public:
*
* Scheme | \f$ x_C \f$ | \f$ t_C \f$ | \f$ \alpha \f$ | \f$ x_N \f$ | \f$ x_O \f$
* -------------- | --------------- | --------------- | --------------------- | --------------- | ----------------------
* Forward Euler | \f$ x_n \f$ | \f$ t_n \f$ | \f$ 1/\Delta t \f$ | \f$ x_{n+1} \f$ | \f$ x_n / \Delta t \f$
* Backward Euler | \f$ x_{n+1} \f$ | \f$ t_{n+1} \f$ | \f$ 1/\Delta t \f$ | \f$ x_{n+1} \f$ | \f$ x_n / \Delta t \f$
* Crank-Nicolson | \f$ x_{n+1} \f$ | \f$ t_{n+1} \f$ | \f$ 1/\Delta t \f$ | \f$ x_{n+1} \f$ | \f$ x_n / \Delta t \f$
* BDF(2) | \f$ x_{n+2} \f$ | \f$ t_{n+1} \f$ | \f$ 3/(2\Delta t) \f$ | \f$ x_{n+2} \f$ | \f$ (2\cdot x_{n+1} - x_n/2)/\Delta t \f$
*
* The other backward differentiation formulas of orders 1 to 6 are also implemented, but only
* BDF(2) has bee given here for brevity.
*/
class TimeDiscretization
{
......@@ -136,11 +109,8 @@ public:
*
* \param t The current timestep.
* \param x The solution at the current timestep.
* \param strg Trigger storing some internal state.
* Currently only used by the CrankNicolson scheme.
*/
virtual void pushState(const double t, GlobalVector const& x,
InternalMatrixStorage const& strg) = 0;
virtual void pushState(const double t, GlobalVector const& x) = 0;
/*!
* Restores the given vector x to its old value.
......@@ -191,43 +161,6 @@ public:
virtual ~TimeDiscretization() = default;
//! \name Extended Interface
//! These methods are provided primarily to make certain concrete time
//! discretizations
//! with special demands, such as the forward Euler or Crank-Nicolson
//! schemes, possible.
//! @{
/*! Tell whether this scheme inherently requires a nonlinear solver or not.
*
* The ForwardEuler scheme is inherently linear in that sense, the others
* are not.
*/
virtual bool isLinearTimeDisc() const { return false; }
/*! Returns \f$ \partial x_C / \partial x_N \f$.
*
* The ForwardEuler scheme overrides this.
*/
virtual double getDxDx() const { return 1.0; }