Skip to content
Snippets Groups Projects
Commit f5b59e9d authored by Norihiro Watanabe's avatar Norihiro Watanabe
Browse files

rm Mat/Vec templates from NonlinearSystem and NonlinearSolver

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