diff --git a/NumLib/ODESolver/NonlinearSolver-impl.h b/NumLib/ODESolver/NonlinearSolver-impl.h
index 16765defdf436dd1816d97cc4bd8c2ee41d1fe2c..7f34a5be5426ba61d6cf53cc9ab2c9796a172ef8 100644
--- a/NumLib/ODESolver/NonlinearSolver-impl.h
+++ b/NumLib/ODESolver/NonlinearSolver-impl.h
@@ -22,26 +22,26 @@
namespace NumLib
{
-template <typename Matrix, typename Vector>
-void NonlinearSolver<Matrix, Vector, NonlinearSolverTag::Picard>::assemble(
- Vector const& x) const
+inline void NonlinearSolver<NonlinearSolverTag::Picard>::assemble(
+ GlobalVector const& x) const
{
_equation_system->assembleMatricesPicard(x);
}
-template <typename Matrix, typename Vector>
-bool NonlinearSolver<Matrix, Vector, NonlinearSolverTag::Picard>::solve(
- Vector& x,
- std::function<void(unsigned, Vector const&)> const& postIterationCallback)
+inline bool NonlinearSolver<NonlinearSolverTag::Picard>::solve(
+ GlobalVector& x,
+ std::function<void(unsigned, GlobalVector const&)> const& postIterationCallback)
{
namespace BLAS = MathLib::BLAS;
auto& sys = *_equation_system;
- auto& A = MathLib::GlobalMatrixProvider<Matrix>::provider.getMatrix(_A_id);
- auto& rhs =
- MathLib::GlobalVectorProvider<Vector>::provider.getVector(_rhs_id);
+ auto& A =
+ MathLib::GlobalMatrixProvider<GlobalMatrix>::provider.getMatrix(_A_id);
+ auto& rhs = MathLib::GlobalVectorProvider<GlobalVector>::provider.getVector(
+ _rhs_id);
auto& x_new =
- MathLib::GlobalVectorProvider<Vector>::provider.getVector(_x_new_id);
+ MathLib::GlobalVectorProvider<GlobalVector>::provider.getVector(
+ _x_new_id);
bool error_norms_met = false;
@@ -136,16 +136,15 @@ bool NonlinearSolver<Matrix, Vector, NonlinearSolverTag::Picard>::solve(
_maxiter);
}
- MathLib::GlobalMatrixProvider<Matrix>::provider.releaseMatrix(A);
- MathLib::GlobalVectorProvider<Vector>::provider.releaseVector(rhs);
- MathLib::GlobalVectorProvider<Vector>::provider.releaseVector(x_new);
+ MathLib::GlobalMatrixProvider<GlobalMatrix>::provider.releaseMatrix(A);
+ MathLib::GlobalVectorProvider<GlobalVector>::provider.releaseVector(rhs);
+ MathLib::GlobalVectorProvider<GlobalVector>::provider.releaseVector(x_new);
return error_norms_met;
}
-template <typename Matrix, typename Vector>
-void NonlinearSolver<Matrix, Vector, NonlinearSolverTag::Newton>::assemble(
- Vector const& x) const
+inline void NonlinearSolver<NonlinearSolverTag::Newton>::assemble(
+ GlobalVector const& x) const
{
_equation_system->assembleResidualNewton(x);
// TODO if the equation system would be reset to nullptr after each
@@ -153,20 +152,20 @@ void NonlinearSolver<Matrix, Vector, NonlinearSolverTag::Newton>::assemble(
// equation every time and could not forget it.
}
-template <typename Matrix, typename Vector>
-bool NonlinearSolver<Matrix, Vector, NonlinearSolverTag::Newton>::solve(
- Vector& x,
- std::function<void(unsigned, Vector const&)> const& postIterationCallback)
+inline bool NonlinearSolver<NonlinearSolverTag::Newton>::solve(
+ GlobalVector& x,
+ std::function<void(unsigned, GlobalVector const&)> const& postIterationCallback)
{
namespace BLAS = MathLib::BLAS;
auto& sys = *_equation_system;
- auto& res =
- MathLib::GlobalVectorProvider<Vector>::provider.getVector(_res_id);
+ auto& res = MathLib::GlobalVectorProvider<GlobalVector>::provider.getVector(
+ _res_id);
auto& minus_delta_x =
- MathLib::GlobalVectorProvider<Vector>::provider.getVector(
+ MathLib::GlobalVectorProvider<GlobalVector>::provider.getVector(
_minus_delta_x_id);
- auto& J = MathLib::GlobalMatrixProvider<Matrix>::provider.getMatrix(_J_id);
+ auto& J =
+ MathLib::GlobalMatrixProvider<GlobalMatrix>::provider.getMatrix(_J_id);
bool error_norms_met = false;
@@ -203,7 +202,7 @@ bool NonlinearSolver<Matrix, Vector, NonlinearSolverTag::Newton>::solve(
// cf.
// http://www.mcs.anl.gov/petsc/petsc-current/docs/manualpages/Vec/VecWAXPY.html
auto& x_new =
- MathLib::GlobalVectorProvider<Vector>::provider.getVector(
+ MathLib::GlobalVectorProvider<GlobalVector>::provider.getVector(
x, _x_new_id);
BLAS::axpy(x_new, -_alpha, minus_delta_x);
@@ -225,7 +224,7 @@ bool NonlinearSolver<Matrix, Vector, NonlinearSolverTag::Newton>::solve(
"iteration"
" has to be repeated.");
// TODO introduce some onDestroy hook.
- MathLib::GlobalVectorProvider<Vector>::provider
+ MathLib::GlobalVectorProvider<GlobalVector>::provider
.releaseVector(x_new);
continue; // That throws the iteration result away.
}
@@ -233,7 +232,7 @@ bool NonlinearSolver<Matrix, Vector, NonlinearSolverTag::Newton>::solve(
// TODO could be done via swap. Note: that also requires swapping
// the ids. Same for the Picard scheme.
BLAS::copy(x_new, x); // copy new solution to x
- MathLib::GlobalVectorProvider<Vector>::provider.releaseVector(
+ MathLib::GlobalVectorProvider<GlobalVector>::provider.releaseVector(
x_new);
}
@@ -272,21 +271,20 @@ bool NonlinearSolver<Matrix, Vector, NonlinearSolverTag::Newton>::solve(
_maxiter);
}
- MathLib::GlobalMatrixProvider<Matrix>::provider.releaseMatrix(J);
- MathLib::GlobalVectorProvider<Vector>::provider.releaseVector(res);
- MathLib::GlobalVectorProvider<Vector>::provider.releaseVector(
+ MathLib::GlobalMatrixProvider<GlobalMatrix>::provider.releaseMatrix(J);
+ MathLib::GlobalVectorProvider<GlobalVector>::provider.releaseVector(res);
+ MathLib::GlobalVectorProvider<GlobalVector>::provider.releaseVector(
minus_delta_x);
return error_norms_met;
}
-template <typename Matrix, typename Vector>
-std::pair<std::unique_ptr<NonlinearSolverBase<Matrix, Vector>>,
- NonlinearSolverTag>
-createNonlinearSolver(MathLib::LinearSolver<Matrix, Vector>& linear_solver,
- BaseLib::ConfigTree const& config)
+inline std::pair<std::unique_ptr<NonlinearSolverBase>, NonlinearSolverTag>
+createNonlinearSolver(
+ MathLib::LinearSolver<GlobalMatrix, GlobalVector>& linear_solver,
+ BaseLib::ConfigTree const& config)
{
- using AbstractNLS = NonlinearSolverBase<Matrix, Vector>;
+ using AbstractNLS = NonlinearSolverBase;
//! \ogs_file_param{prj__nonlinear_solvers__nonlinear_solver__type}
auto const type = config.getConfigParameter<std::string>("type");
@@ -298,7 +296,7 @@ createNonlinearSolver(MathLib::LinearSolver<Matrix, Vector>& linear_solver,
if (type == "Picard")
{
auto const tag = NonlinearSolverTag::Picard;
- using ConcreteNLS = NonlinearSolver<Matrix, Vector, tag>;
+ using ConcreteNLS = NonlinearSolver<tag>;
return std::make_pair(std::unique_ptr<AbstractNLS>(new ConcreteNLS{
linear_solver, tol, max_iter}),
tag);
@@ -306,7 +304,7 @@ createNonlinearSolver(MathLib::LinearSolver<Matrix, Vector>& linear_solver,
else if (type == "Newton")
{
auto const tag = NonlinearSolverTag::Newton;
- using ConcreteNLS = NonlinearSolver<Matrix, Vector, tag>;
+ using ConcreteNLS = NonlinearSolver<tag>;
return std::make_pair(std::unique_ptr<AbstractNLS>(new ConcreteNLS{
linear_solver, tol, max_iter}),
tag);
diff --git a/NumLib/ODESolver/NonlinearSolver.h b/NumLib/ODESolver/NonlinearSolver.h
index 4ce05ba5c53473f5bfa2215dd1450db2dad9a8f9..140fe406c230c63dcbfee3c9941b8b2348327253 100644
--- a/NumLib/ODESolver/NonlinearSolver.h
+++ b/NumLib/ODESolver/NonlinearSolver.h
@@ -36,7 +36,6 @@ namespace NumLib
* equation.
* \tparam Vector the type of the solution vector of the equation.
*/
-template <typename Matrix, typename Vector>
class NonlinearSolverBase
{
public:
@@ -48,7 +47,7 @@ public:
*
* \param x the state at which the equation system will be assembled.
*/
- virtual void assemble(Vector const& x) const = 0;
+ virtual void assemble(GlobalVector const& x) const = 0;
/*! Assemble and solve the equation system.
*
@@ -58,8 +57,8 @@ public:
* \retval true if the equation system could be solved
* \retval false otherwise
*/
- virtual bool solve(Vector& x,
- std::function<void(unsigned, Vector const&)> const&
+ virtual bool solve(GlobalVector& x,
+ std::function<void(unsigned, GlobalVector const&)> const&
postIterationCallback) = 0;
virtual ~NonlinearSolverBase() = default;
@@ -75,7 +74,7 @@ public:
* \tparam Vector the type of the solution vector of the equation.
* \tparam NLTag a tag indicating the method used for solving the equation.
*/
-template <typename Matrix, typename Vector, NonlinearSolverTag NLTag>
+template <NonlinearSolverTag NLTag>
class NonlinearSolver;
/*! Find a solution to a nonlinear equation using the Newton-Raphson method.
@@ -84,14 +83,14 @@ class NonlinearSolver;
* equation.
* \tparam Vector the type of the solution vector of the equation.
*/
-template <typename Matrix, typename Vector>
-class NonlinearSolver<Matrix, Vector, NonlinearSolverTag::Newton> final
- : public NonlinearSolverBase<Matrix, Vector>
+template <>
+class NonlinearSolver<NonlinearSolverTag::Newton> final
+ : public NonlinearSolverBase
{
public:
//! Type of the nonlinear equation system to be solved.
- using System = NonlinearSystem<Matrix, Vector, NonlinearSolverTag::Newton>;
- using LinearSolver = MathLib::LinearSolver<Matrix, Vector>;
+ using System = NonlinearSystem<NonlinearSolverTag::Newton>;
+ using LinearSolver = MathLib::LinearSolver<GlobalMatrix, GlobalVector>;
/*! Constructs a new instance.
*
@@ -109,10 +108,10 @@ public:
//! Set the nonlinear equation system that will be solved.
void setEquationSystem(System& eq) { _equation_system = &eq; }
- void assemble(Vector const& x) const override;
+ void assemble(GlobalVector const& x) const override;
- bool solve(Vector& x,
- std::function<void(unsigned, Vector const&)> const&
+ bool solve(GlobalVector& x,
+ std::function<void(unsigned, GlobalVector const&)> const&
postIterationCallback) override;
private:
@@ -139,14 +138,14 @@ private:
* equation.
* \tparam Vector the type of the solution vector of the equation.
*/
-template <typename Matrix, typename Vector>
-class NonlinearSolver<Matrix, Vector, NonlinearSolverTag::Picard> final
- : public NonlinearSolverBase<Matrix, Vector>
+template <>
+class NonlinearSolver<NonlinearSolverTag::Picard> final
+ : public NonlinearSolverBase
{
public:
//! Type of the nonlinear equation system to be solved.
- using System = NonlinearSystem<Matrix, Vector, NonlinearSolverTag::Picard>;
- using LinearSolver = MathLib::LinearSolver<Matrix, Vector>;
+ using System = NonlinearSystem<NonlinearSolverTag::Picard>;
+ using LinearSolver = MathLib::LinearSolver<GlobalMatrix, GlobalVector>;
/*! Constructs a new instance.
*
@@ -164,10 +163,10 @@ public:
//! Set the nonlinear equation system that will be solved.
void setEquationSystem(System& eq) { _equation_system = &eq; }
- void assemble(Vector const& x) const override;
+ void assemble(GlobalVector const& x) const override;
- bool solve(Vector& x,
- std::function<void(unsigned, Vector const&)> const&
+ bool solve(GlobalVector& x,
+ std::function<void(unsigned, GlobalVector const&)> const&
postIterationCallback) override;
private:
@@ -193,11 +192,10 @@ private:
* nonlinear solver instance and the \c tag indicates if it uses
* the Picard or Newton-Raphson method
*/
-template <typename Matrix, typename Vector>
-std::pair<std::unique_ptr<NonlinearSolverBase<Matrix, Vector>>,
- NonlinearSolverTag>
-createNonlinearSolver(MathLib::LinearSolver<Matrix, Vector>& linear_solver,
- BaseLib::ConfigTree const& config);
+std::pair<std::unique_ptr<NonlinearSolverBase>, NonlinearSolverTag>
+createNonlinearSolver(
+ MathLib::LinearSolver<GlobalMatrix, GlobalVector>& linear_solver,
+ BaseLib::ConfigTree const& config);
//! @}
diff --git a/NumLib/ODESolver/NonlinearSystem.h b/NumLib/ODESolver/NonlinearSystem.h
index ecbaef309615932116c6c1072bd844afa0090190..c1eff192837733bc7205b4116b919819b01a9ce3 100644
--- a/NumLib/ODESolver/NonlinearSystem.h
+++ b/NumLib/ODESolver/NonlinearSystem.h
@@ -25,7 +25,7 @@ namespace NumLib
* \tparam Vector the type of the solution vector of the equation.
* \tparam NLTag a tag indicating the method used for solving the equation.
*/
-template <typename Matrix, typename Vector, NonlinearSolverTag NLTag>
+template <NonlinearSolverTag NLTag>
class NonlinearSystem;
/*! A System of nonlinear equations to be solved with the Newton-Raphson method.
@@ -37,16 +37,15 @@ class NonlinearSystem;
* equation.
* \tparam Vector the type of the solution vector of the equation.
*/
-template <typename Matrix, typename Vector>
-class NonlinearSystem<Matrix, Vector, NonlinearSolverTag::Newton>
- : public EquationSystem<Vector>
+template <>
+class NonlinearSystem<NonlinearSolverTag::Newton> : public EquationSystem
{
public:
//! Assembles the residual at the point \c x.
- virtual void assembleResidualNewton(Vector const& x) = 0;
+ virtual void assembleResidualNewton(GlobalVector const& x) = 0;
//! Assembles the Jacobian of the residual at the point \c x.
- virtual void assembleJacobian(Vector const& x) = 0;
+ virtual void assembleJacobian(GlobalVector const& x) = 0;
/*! Writes the residual at point \c x to \c res.
*
@@ -55,21 +54,22 @@ public:
*
* \todo Remove argument \c x.
*/
- virtual void getResidual(Vector const& x, Vector& res) const = 0;
+ virtual void getResidual(GlobalVector const& x,
+ GlobalVector& res) const = 0;
/*! Writes the Jacobian of the residual to \c Jac.
*
* \pre assembleJacobian() must have been called before.
*/
- virtual void getJacobian(Matrix& Jac) const = 0;
+ virtual void getJacobian(GlobalMatrix& Jac) const = 0;
//! Apply known solutions to the solution vector \c x.
- virtual void applyKnownSolutions(Vector& x) const = 0;
+ virtual void applyKnownSolutions(GlobalVector& x) const = 0;
//! Apply known solutions to the linearized equation system
//! \f$ \mathit{Jac} \cdot (-\Delta x) = \mathit{res} \f$.
- virtual void applyKnownSolutionsNewton(Matrix& Jac, Vector& res,
- Vector& minus_delta_x) = 0;
+ virtual void applyKnownSolutionsNewton(GlobalMatrix& Jac, GlobalVector& res,
+ GlobalVector& minus_delta_x) = 0;
};
/*! A System of nonlinear equations to be solved with the Picard fixpoint
@@ -82,27 +82,26 @@ public:
* equation.
* \tparam Vector the type of the solution vector of the equation.
*/
-template <typename Matrix, typename Vector>
-class NonlinearSystem<Matrix, Vector, NonlinearSolverTag::Picard>
- : public EquationSystem<Vector>
+template <>
+class NonlinearSystem<NonlinearSolverTag::Picard> : public EquationSystem
{
public:
//! Assembles the linearized equation at point \c x.
- virtual void assembleMatricesPicard(Vector const& x) = 0;
+ virtual void assembleMatricesPicard(GlobalVector const& x) = 0;
//! Writes the linearized equation system matrix to \c A.
- virtual void getA(Matrix& A) const = 0;
+ virtual void getA(GlobalMatrix& A) const = 0;
//! Writes the linearized equation system right-hand side to \c rhs.
- virtual void getRhs(Vector& rhs) const = 0;
+ virtual void getRhs(GlobalVector& rhs) const = 0;
//! Apply known solutions to the solution vector \c x.
- virtual void applyKnownSolutions(Vector& x) const = 0;
+ virtual void applyKnownSolutions(GlobalVector& x) const = 0;
//! Apply known solutions to the linearized equation system
//! \f$ A \cdot x = \mathit{rhs} \f$.
- virtual void applyKnownSolutionsPicard(Matrix& A, Vector& rhs,
- Vector& x) = 0;
+ virtual void applyKnownSolutionsPicard(GlobalMatrix& A, GlobalVector& rhs,
+ GlobalVector& x) = 0;
};
//! @}