From 6fafe8bc8541ffeaebfadf1b2a0567c4725f2d1a Mon Sep 17 00:00:00 2001
From: Christoph Lehmann <christoph.lehmann@ufz.de>
Date: Wed, 17 Feb 2016 17:32:31 +0100
Subject: [PATCH] [NL] comments from code review
---
NumLib/ODESolver/MatrixTranslator.h | 1 +
NumLib/ODESolver/NonlinearSolver-impl.h | 17 +++++++++++-----
NumLib/ODESolver/NonlinearSolver.h | 9 +++++----
NumLib/ODESolver/NonlinearSystem.h | 8 +++++++-
NumLib/ODESolver/ODESystem.h | 5 +++--
NumLib/ODESolver/TimeDiscretization.h | 22 +++++++++++++--------
NumLib/ODESolver/TimeDiscretizedODESystem.h | 6 ------
NumLib/ODESolver/TimeLoop.h | 3 +++
8 files changed, 45 insertions(+), 26 deletions(-)
diff --git a/NumLib/ODESolver/MatrixTranslator.h b/NumLib/ODESolver/MatrixTranslator.h
index dc0162c6e20..8ab72fbe636 100644
--- a/NumLib/ODESolver/MatrixTranslator.h
+++ b/NumLib/ODESolver/MatrixTranslator.h
@@ -64,6 +64,7 @@ public:
//! Computes the Jacobian of the residual and writes it to \c Jac_out.
virtual void getJacobian(Matrix const& Jac_in, Matrix& Jac_out) const = 0;
+ // TODO check (void)
/*! Allows to store the given matrices internally for later use.
*
* \remark
diff --git a/NumLib/ODESolver/NonlinearSolver-impl.h b/NumLib/ODESolver/NonlinearSolver-impl.h
index 5fc9b8ed11a..8cf2e92032e 100644
--- a/NumLib/ODESolver/NonlinearSolver-impl.h
+++ b/NumLib/ODESolver/NonlinearSolver-impl.h
@@ -43,12 +43,15 @@ solve(NonlinearSystem<Matrix, Vector, NonlinearSolverTag::Picard> &sys, Vector &
bool success = false;
+ BLAS::copy(x, _x_new); // set initial guess, TODO save the copy
+
for (unsigned iteration=1; iteration<_maxiter; ++iteration)
{
- sys.assembleMatricesPicard(x);
+ sys.assembleMatricesPicard(_x_new);
sys.getA(_A);
sys.getRhs(_rhs);
- // TODO _x_new might not have been initialized
+
+ // Here _x_new has to be used and it has to be equal to x!
sys.applyKnownComponentsPicard(_A, _rhs, _x_new);
// std::cout << "A:\n" << Eigen::MatrixXd(A) << "\n";
@@ -56,12 +59,11 @@ solve(NonlinearSystem<Matrix, Vector, NonlinearSolverTag::Picard> &sys, Vector &
oneShotLinearSolve(_A, _rhs, _x_new);
+ // x is used as delta_x in order to compute the error.
BLAS::aypx(x, -1.0, _x_new); // x = _x_new - x
auto const error = norm(x);
// INFO(" picard iteration %u error: %e", iteration, error);
- x = _x_new;
-
if (error < _tol) {
success = true;
break;
@@ -74,6 +76,8 @@ solve(NonlinearSystem<Matrix, Vector, NonlinearSolverTag::Picard> &sys, Vector &
}
}
+ x = _x_new; // always set to what we computed last
+
return success;
}
@@ -96,6 +100,8 @@ solve(NonlinearSystem<Matrix, Vector, NonlinearSolverTag::Newton> &sys, Vector &
bool success = false;
+ // TODO init _minus_delta_x to right size and 0.0
+
for (unsigned iteration=1; iteration<_maxiter; ++iteration)
{
sys.assembleResidualNewton(x);
@@ -103,6 +109,7 @@ solve(NonlinearSystem<Matrix, Vector, NonlinearSolverTag::Newton> &sys, Vector &
// std::cout << " res:\n" << res << std::endl;
+ // TODO streamline the, make consistent with Picard.
if (norm(_res) < _tol) {
success = true;
break;
@@ -120,7 +127,7 @@ solve(NonlinearSystem<Matrix, Vector, NonlinearSolverTag::Newton> &sys, Vector &
// auto const dx_norm = _minus_delta_x.norm();
// INFO(" newton iteration %u, norm of delta x: %e", iteration, dx_norm);
- BLAS::axpy(x, -1.0, _minus_delta_x);
+ BLAS::axpy(x, -_alpha, _minus_delta_x);
if (sys.isLinear()) {
// INFO(" newton linear system. not looping");
diff --git a/NumLib/ODESolver/NonlinearSolver.h b/NumLib/ODESolver/NonlinearSolver.h
index 02f1ca009d9..896e8d44107 100644
--- a/NumLib/ODESolver/NonlinearSolver.h
+++ b/NumLib/ODESolver/NonlinearSolver.h
@@ -80,10 +80,11 @@ private:
const double _tol; //!< tolerance of the solver
const unsigned _maxiter; //!< maximum number of iterations
- Vector _res; //!< The residual.
- Matrix _J; //!< The Jacobian of the residual.
- Vector _minus_delta_x; //!< The Newton-Raphson method solves the linearized equation
- //!< \f$ J \cdot (-\Delta x) = r \f$ repeatedly.
+ Vector _res; //!< The residual.
+ Matrix _J; //!< The Jacobian of the residual.
+ Vector _minus_delta_x; //!< The Newton-Raphson method solves the linearized equation
+ //!< \f$ J \cdot (-\Delta x) = r \f$ repeatedly.
+ double const _alpha = 1; //!< damping factor
};
diff --git a/NumLib/ODESolver/NonlinearSystem.h b/NumLib/ODESolver/NonlinearSystem.h
index 5af8a865616..c61847442d3 100644
--- a/NumLib/ODESolver/NonlinearSystem.h
+++ b/NumLib/ODESolver/NonlinearSystem.h
@@ -62,6 +62,7 @@ public:
*/
virtual void getJacobian(Matrix& Jac) const = 0;
+ // TODO solutions
//! Apply known solutions to the linearized equation system
//! \f$ \mathit{Jac} \cdot (-\Delta x) = \mathit{res} \f$.
virtual void applyKnownComponentsNewton(
@@ -76,10 +77,12 @@ public:
*/
virtual bool isLinear() const = 0;
+ // TODO add getNumEquations()
+
virtual ~NonlinearSystem() = default;
};
-
+// TODO common base class
/*! A System of nonlinear equations to be solved with the Picard fixpoint
* iteration method.
*
@@ -102,6 +105,7 @@ public:
//! Writes the linearized equation system right-hand side to \c rhs.
virtual void getRhs(Vector& rhs) const = 0;
+ // TODO components
//! Apply known solutions to the linearized equation system
//! \f$ A \cdot x = \mathit{rhs} \f$.
virtual void applyKnownComponentsPicard(
@@ -116,6 +120,8 @@ public:
*/
virtual bool isLinear() const = 0;
+ // TODO add getNumEquations()
+
virtual ~NonlinearSystem() = default;
};
diff --git a/NumLib/ODESolver/ODESystem.h b/NumLib/ODESolver/ODESystem.h
index 37db8a9d30e..9c457754006 100644
--- a/NumLib/ODESolver/ODESystem.h
+++ b/NumLib/ODESolver/ODESystem.h
@@ -73,7 +73,7 @@ public:
using Index = typename MatrixTraits<Matrix>::Index;
- // TODO doc
+ // TODO solutions
virtual std::vector<ProcessLib::DirichletBc<Index> > const* getKnownComponents() const
{
return nullptr; // by default there are no known components
@@ -85,7 +85,8 @@ public:
/*! Interface for a first-order implicit quasi-linear ODE.
*
- * ODEs using this interface also provide a Jacobian.
+ * ODEs using this interface also provide a Jacobian in addition
+ * to the functionality of the Picard-related interface.
*/
template<typename Matrix, typename Vector>
class ODESystem<Matrix, Vector,
diff --git a/NumLib/ODESolver/TimeDiscretization.h b/NumLib/ODESolver/TimeDiscretization.h
index 45efa40c52b..05721c62aa8 100644
--- a/NumLib/ODESolver/TimeDiscretization.h
+++ b/NumLib/ODESolver/TimeDiscretization.h
@@ -153,17 +153,18 @@ public:
auto const dxdot_dx = getCurrentXWeight();
+ // xdot = dxdot_dx * x_at_new_timestep - x_old
getWeightedOldX(xdot);
- BLAS::axpby(xdot, dxdot_dx, -1.0, x_at_new_timestep);
+ BLAS::axpby(xdot, dxdot_dx, -1.0, x_at_new_timestep); // TODO consistent
}
//! Returns \f$ \alpha = \partial \hat x / \partial x_C \f$.
- virtual double getCurrentXWeight() const = 0;
+ virtual double getCurrentXWeight() const = 0; // TODO maybe getNewXWeight
//! Returns \f$ x_O \f$.
virtual void getWeightedOldX(Vector& y) const = 0; // = x_old
- ~TimeDiscretization() = default;
+ virtual ~TimeDiscretization() = default; // TODO virtual needed?
//! \name Extended Interface
//! These methods are provided primarily to make certain concrete time discretizations
@@ -212,9 +213,8 @@ public:
_x_old = x0;
}
- void pushState(const double t, Vector const& x, InternalMatrixStorage const&) override
+ void pushState(const double /*t*/, Vector const& x, InternalMatrixStorage const&) override
{
- (void) t;
_x_old = x;
}
@@ -234,6 +234,8 @@ public:
void getWeightedOldX(Vector& y) const override
{
namespace BLAS = MathLib::BLAS;
+
+ // y = x_old / delta_t
BLAS::copy(_x_old, y);
BLAS::scale(y, 1.0/_delta_t);
}
@@ -256,9 +258,8 @@ public:
_x_old = x0;
}
- void pushState(const double t, Vector const& x, InternalMatrixStorage const&) override
+ void pushState(const double /*t*/, Vector const& x, InternalMatrixStorage const&) override
{
- (void) t;
_x_old = x;
}
@@ -283,6 +284,8 @@ public:
void getWeightedOldX(Vector& y) const override
{
namespace BLAS = MathLib::BLAS;
+
+ // y = x_old / delta_t
BLAS::copy(_x_old, y);
BLAS::scale(y, 1.0/_delta_t);
}
@@ -351,6 +354,8 @@ public:
void getWeightedOldX(Vector& y) const override
{
namespace BLAS = MathLib::BLAS;
+
+ // y = x_old / delta_t
BLAS::copy(_x_old, y);
BLAS::scale(y, 1.0/_delta_t);
}
@@ -384,7 +389,7 @@ const double BDF_Coeffs[6][7] = {
{ 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, 0.2 },
+ { 137.0 / 60.0, 5.0, -5.0, 10.0 / 3.0, -1.25, 0.2 }, // TODO change to rational numbers
{ 147.0 / 60.0, 6.0, -7.5, 20.0 / 3.0, -3.75, 1.2, -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
@@ -424,6 +429,7 @@ public:
void pushState(const double t, Vector const& x, InternalMatrixStorage const&) override
{
(void) t;
+ // TODO use boost cirular buffer?
// until _xs_old is filled, lower-order BDF formulas are used.
if (_xs_old.size() < _num_steps) {
diff --git a/NumLib/ODESolver/TimeDiscretizedODESystem.h b/NumLib/ODESolver/TimeDiscretizedODESystem.h
index 42be649191c..ef7f70b8653 100644
--- a/NumLib/ODESolver/TimeDiscretizedODESystem.h
+++ b/NumLib/ODESolver/TimeDiscretizedODESystem.h
@@ -133,12 +133,6 @@ public:
void getResidual(Vector const& x_new_timestep, Vector& res) const override
{
- // TODO Maybe the duplicate calculation of xdot here and in assembleJacobian
- // can be optimuized. However, that would make the interface a bit more
- // fragile.
- _time_disc.getXdot(x_new_timestep, _xdot);
-
- // TODO change to calculate
_mat_trans->getResidual(_M, _K, _b, x_new_timestep, _xdot, res);
}
diff --git a/NumLib/ODESolver/TimeLoop.h b/NumLib/ODESolver/TimeLoop.h
index 699c9d26796..60d0fad9447 100644
--- a/NumLib/ODESolver/TimeLoop.h
+++ b/NumLib/ODESolver/TimeLoop.h
@@ -68,6 +68,7 @@ private:
//! @}
+// TODO rename class
template<typename Matrix, typename Vector, NonlinearSolverTag NLTag>
template<typename Callback>
bool
@@ -86,6 +87,7 @@ loop(const double t0, const Vector x0, const double t_end, const double delta_t,
time_disc.pushState(t0, x0, _ode_sys); // TODO: that might do duplicate work
}
+ // TODO change n*delta_t
double t;
unsigned timestep = 0;
bool nl_slv_succeeded = true;
@@ -99,6 +101,7 @@ loop(const double t0, const Vector x0, const double t_end, const double delta_t,
time_disc.pushState(t, x, _ode_sys);
+ // TODO check const ref
auto const t_cb = t; // make sure the callback cannot overwrite anything.
auto const& x_cb = x; // ditto.
post_timestep(t_cb, x_cb);
--
GitLab