https://gitlab.opengeosys.org/ogs/ogs.git
Tip revision: 0fde67c0ca8818233131619e1c07b8c4be4094a0 authored by Christoph Lehmann on 01 May 2016, 20:11:48 UTC
[TES] chaged reaction taming strategy
[TES] chaged reaction taming strategy
Tip revision: 0fde67c
ConcreteODESolver.h
/**
* \copyright
* Copyright (c) 2012-2016, OpenGeoSys Community (http://www.opengeosys.org)
* Distributed under a Modified BSD License.
* See accompanying file LICENSE.txt or
* http://www.opengeosys.org/project/license
*
*/
#ifndef MATHLIB_ODE_CONCRETEODESOLVER_H
#define MATHLIB_ODE_CONCRETEODESOLVER_H
#include <memory>
#include "FunctionHandles.h"
#include "ODESolver.h"
namespace BaseLib { class ConfigTree; }
namespace MathLib
{
namespace ODE
{
template <unsigned NumEquations>
std::unique_ptr<ODESolver<NumEquations>> createODESolver(
BaseLib::ConfigTree const& config);
//! \addtogroup ExternalODESolverInterface
//! @{
/*! ODE solver with a bounds-safe interface using a concrete implementation.
*
* This class makes contact between the abstract \c ODESolver interface and a
* certain solver \c Implementation.
*
* The interface of this class inherits the array bounds checking from \c
* ODESolver.
* Its methods forward calls to the \c Implementation erasing array bounds info
* by passing \c std::array and Eigen::Matrix arguments as raw pointers.
*
* Thus the \c Implementation does not need template to be templated
* which makes it possible for \c Implementation to use the pimpl idiom whereby
* the headers of external libraries only have to be included in the final cpp
* file. Thereby our namespaces do not get polluted with symbols from external
* libraries.
*
* \tparam NumEquations the number of equations in the ODE system.
* \tparam Implementation a concrete ODE solver implementation used as a backend.
*/
template <typename Implementation,
unsigned NumEquations>
class ConcreteODESolver final
: public ODESolver<NumEquations>,
private Implementation
{
public:
void setFunction(Function<NumEquations> f,
JacobianFunction<NumEquations> df) override
{
Implementation::setFunction(
std::unique_ptr<
detail::FunctionHandlesImpl<NumEquations>>{
new detail::FunctionHandlesImpl<NumEquations>{f, df}});
}
void setTolerance(const std::array<double, NumEquations>& abstol,
const double reltol) override
{
Implementation::setTolerance(abstol.data(), reltol);
}
void setTolerance(const double abstol, const double reltol) override
{
Implementation::setTolerance(abstol, reltol);
}
void setIC(const double t0,
std::initializer_list<double> const& y0) override
{
assert(y0.size() == NumEquations);
Implementation::setIC(t0, y0.begin());
}
void setIC(const double t0,
Eigen::Matrix<double, NumEquations, 1, Eigen::ColMajor> const& y0
) override
{
Implementation::setIC(t0, y0.data());
}
void preSolve() override { Implementation::preSolve(); }
bool solve(const double t) override { return Implementation::solve(t); }
MappedConstVector<NumEquations> getSolution() const override
{
return MappedConstVector<NumEquations>{Implementation::getSolution()};
}
double getTime() const override { return Implementation::getTime(); }
Eigen::Matrix<double, NumEquations, 1, Eigen::ColMajor> getYDot(
const double t, const MappedConstVector<NumEquations>& y) const override
{
Eigen::Matrix<double, NumEquations, 1, Eigen::ColMajor> y_dot;
Implementation::getYDot(t, y.data(), y_dot.data());
return y_dot;
}
private:
//! Instances of this class shall only be constructed by createODESolver().
ConcreteODESolver(BaseLib::ConfigTree const& config)
: Implementation{config, NumEquations}
{
}
friend std::unique_ptr<ODESolver<NumEquations>>
createODESolver<NumEquations>(
BaseLib::ConfigTree const& config);
};
//! @}
} // namespace ODE
} // namespace MathLib
#endif // MATHLIB_ODE_CONCRETEODESOLVER_H