From 49ba9606a7cdcf74e2f1a75613412beeaa3bfbeb Mon Sep 17 00:00:00 2001
From: Wenqing Wang <wenqing.wang@ufz.de>
Date: Fri, 28 Apr 2017 17:41:58 +0200
Subject: [PATCH] [dt] Introduced a member _dx to store the solution increment

as a buffer for computing the relative error.
---
 NumLib/ODESolver/TimeDiscretization.cpp | 30 ++++++++++++-------------
 NumLib/ODESolver/TimeDiscretization.h   | 29 ++++++++++++++++++++----
 2 files changed, 40 insertions(+), 19 deletions(-)

diff --git a/NumLib/ODESolver/TimeDiscretization.cpp b/NumLib/ODESolver/TimeDiscretization.cpp
index a3236ca1c5c..f25f6371215 100644
--- a/NumLib/ODESolver/TimeDiscretization.cpp
+++ b/NumLib/ODESolver/TimeDiscretization.cpp
@@ -10,20 +10,14 @@
  */
 
 #include "TimeDiscretization.h"
+
+#include "MathLib/LinAlg/MatrixVectorTraits.h"
+
 namespace NumLib
 {
-/**
- * Compute and return the relative change of solutions between two successive
- * time steps by \f$ e_n = \|u^{n+1}-u^{n}\|/\|u^{n+1}\| \f$.
- *
- * @param x      The current solution
- * @param x_old  The previous solution
- * @return       \f$ e_n = \|u^{n+1}-u^{n}\|/\|u^{n+1}\| \f$.
- *
- *  \warning the value of x_old is changed to x - x_old after this computation.
- */
-static double computeRelativeError(GlobalVector const& x, GlobalVector& x_old,
-                                   MathLib::VecNormType norm_type)
+double TimeDiscretization::computeRelativeError(GlobalVector const& x,
+                                                GlobalVector const& x_old,
+                                                MathLib::VecNormType norm_type)
 {
     if (norm_type == MathLib::VecNormType::INVALID)
         return 0.;
@@ -31,9 +25,15 @@ static double computeRelativeError(GlobalVector const& x, GlobalVector& x_old,
     const double norm2_x = MathLib::LinAlg::norm(x, norm_type);
     assert(norm2_x > std::numeric_limits<double>::epsilon());
 
-    // dx = x - x_old --> x_old
-    MathLib::LinAlg::axpy(x_old, -1.0, x);
-    return MathLib::LinAlg::norm(x_old, norm_type) / norm2_x;
+    if (!_dx)
+        _dx = MathLib::MatrixVectorTraits<GlobalVector>::newInstance(x);
+   
+    auto& dx = *_dx;
+    MathLib::LinAlg::copy(x_old, dx); // copy x_old to dx.
+
+    // dx = x - x_old --> x - dx --> dx
+    MathLib::LinAlg::axpy(dx, -1.0, x);
+    return MathLib::LinAlg::norm(dx, norm_type) / norm2_x;
 }
 
 double BackwardEuler::getRelativeError(GlobalVector const& x,
diff --git a/NumLib/ODESolver/TimeDiscretization.h b/NumLib/ODESolver/TimeDiscretization.h
index 6f2044d6c9f..930b5b4f876 100644
--- a/NumLib/ODESolver/TimeDiscretization.h
+++ b/NumLib/ODESolver/TimeDiscretization.h
@@ -114,6 +114,8 @@ public:
 class TimeDiscretization
 {
 public:
+    TimeDiscretization() : _dx(nullptr) {}
+
     //! Sets the initial condition.
     virtual void setInitialState(const double t0, GlobalVector const& x0) = 0;
 
@@ -211,6 +213,23 @@ public:
      */
     virtual bool needsPreload() const { return false; }
     //! @}
+
+protected:
+    std::unique_ptr<GlobalVector> _dx; ///< Used to store \f$ u_{n+1}-u_{n}\f$.
+
+    /**
+      * Compute and return the relative change of solutions between two successive
+      * time steps by \f$ e_n = \|u^{n+1}-u^{n}\|/\|u^{n+1}\| \f$.
+      *
+      * @param x      The current solution
+      * @param x_old  The previous solution
+      * @return       \f$ e_n = \|u^{n+1}-u^{n}\|/\|u^{n+1}\| \f$.
+      *
+      *  \warning the value of x_old is changed to x - x_old after this computation.
+    */    
+    double computeRelativeError(GlobalVector const& x,
+                                GlobalVector const& x_old,
+                                MathLib::VecNormType norm_type);    
 };
 
 //! Backward Euler scheme.
@@ -218,7 +237,8 @@ class BackwardEuler final : public TimeDiscretization
 {
 public:
     BackwardEuler()
-        : _x_old(NumLib::GlobalVectorProvider::provider.getVector())
+        : TimeDiscretization(),
+          _x_old(NumLib::GlobalVectorProvider::provider.getVector())
     {
     }
 
@@ -270,7 +290,8 @@ class ForwardEuler final : public TimeDiscretization
 {
 public:
     ForwardEuler()
-        : _x_old(NumLib::GlobalVectorProvider::provider.getVector())
+        :  TimeDiscretization(),
+           _x_old(NumLib::GlobalVectorProvider::provider.getVector())
     {
     }
 
@@ -347,7 +368,7 @@ public:
      *              \arg 0.5 traditional Crank-Nicolson scheme.
      */
     explicit CrankNicolson(const double theta)
-        : _theta(theta),
+        : TimeDiscretization(), _theta(theta),
           _x_old(NumLib::GlobalVectorProvider::provider.getVector())
     {
     }
@@ -419,7 +440,7 @@ public:
      *       the first timesteps.
      */
     explicit BackwardDifferentiationFormula(const unsigned num_steps)
-        : _num_steps(num_steps)
+        :  TimeDiscretization(), _num_steps(num_steps)
     {
         assert(1 <= num_steps && num_steps <= 6);
         _xs_old.reserve(num_steps);
-- 
GitLab