OGS 6.2.0-97-g4a610c866
TimeDiscretization.cpp
Go to the documentation of this file.
1 
12 #include "TimeDiscretization.h"
13 
15 
16 namespace NumLib
17 {
19  GlobalVector const& x,
20  GlobalVector const& x_old,
21  MathLib::VecNormType norm_type)
22 {
23  if (norm_type == MathLib::VecNormType::INVALID)
24  {
25  OGS_FATAL("An invalid norm type has been passed");
26  }
27 
28  if (!_dx)
29  {
31  }
32 
33  auto& dx = *_dx;
34  MathLib::LinAlg::copy(x, dx); // copy x to dx.
35 
36  // dx = x - x_old --> x - dx --> dx
37  MathLib::LinAlg::axpy(dx, -1.0, x_old);
38  const double norm_dx = MathLib::LinAlg::norm(dx, norm_type);
39 
40  const double norm_x = MathLib::LinAlg::norm(x, norm_type);
41  if (norm_x > std::numeric_limits<double>::epsilon())
42  {
43  return norm_dx / norm_x;
44  }
45 
46  // Both of norm_x and norm_dx are close to zero
47  if (norm_dx < std::numeric_limits<double>::epsilon())
48  {
49  return 1.0;
50  }
51 
52  // Only norm_x is close to zero
53  return norm_dx / std::numeric_limits<double>::epsilon();
54 }
55 
57  GlobalVector const& x, MathLib::VecNormType norm_type)
58 {
59  return computeRelativeChangeFromPreviousTimestep(x, _x_old, norm_type);
60 }
61 
63  GlobalVector const& x, MathLib::VecNormType norm_type)
64 {
65  return computeRelativeChangeFromPreviousTimestep(x, _x_old, norm_type);
66 }
67 
69  GlobalVector const& x, MathLib::VecNormType norm_type)
70 {
71  return computeRelativeChangeFromPreviousTimestep(x, _x_old, norm_type);
72 }
73 
75  GlobalVector const& x, MathLib::VecNormType norm_type)
76 {
78  x, *_xs_old[_offset], norm_type);
79 }
80 
82  GlobalVector const& x,
83  InternalMatrixStorage const&)
84 {
85  namespace LinAlg = MathLib::LinAlg;
86  // TODO use boost circular buffer?
87 
88  // until _xs_old is filled, lower-order BDF formulas are used.
89  if (_xs_old.size() < _num_steps)
90  {
91  _xs_old.push_back(&NumLib::GlobalVectorProvider::provider.getVector(x));
92  }
93  else
94  {
95  LinAlg::copy(x, *_xs_old[_offset]);
96  _offset =
97  (_offset + 1) % _num_steps; // treat _xs_old as a circular buffer
98  }
99 }
100 
101 namespace detail
102 {
104 const double BDF_Coeffs[6][7] = {
105  // leftmost column: weight of the solution at the new timestep
106  // signs of columns > 1 are flipped compared to standard BDF tableaus
107  {1.0, 1.0},
108  {1.5, 2.0, -0.5},
109  {11.0 / 6.0, 3.0, -1.5, 1.0 / 3.0},
110  {25.0 / 12.0, 4.0, -3.0, 4.0 / 3.0, -0.25},
111  {137.0 / 60.0, 5.0, -5.0, 10.0 / 3.0, -1.25, 1.0 / 5.0},
112  {147.0 / 60.0, 6.0, -7.5, 20.0 / 3.0, -3.75, 6.0 / 5.0, -1.0 / 6.0}
113  // coefficient of (for BDF(6), the oldest state, x_n, is always rightmost)
114  // x_+6, x_+5, x_+4, x_+3, x_+2, x_+1, x_n
115 };
116 } // namespace detail
117 
119 {
120  auto const k = eff_num_steps();
121  return detail::BDF_Coeffs[k - 1][0] / _delta_t;
122 }
123 
125 {
126  namespace LinAlg = MathLib::LinAlg;
127 
128  auto const k = eff_num_steps();
129  auto const* const BDFk = detail::BDF_Coeffs[k - 1];
130 
131  // compute linear combination \sum_{i=0}^{k-1} BDFk_{k-i} \cdot x_{n+i}
132  LinAlg::copy(*_xs_old[_offset], y); // _xs_old[offset] = x_n
133  LinAlg::scale(y, BDFk[k]);
134 
135  for (unsigned i = 1; i < k; ++i)
136  {
137  auto const off = (_offset + i) % k;
138  LinAlg::axpy(y, BDFk[k - i], *_xs_old[off]);
139  }
140 
141  LinAlg::scale(y, 1.0 / _delta_t);
142 }
143 
144 } // end of namespace NumLib
double getRelativeChangeFromPreviousTimestep(GlobalVector const &x, MathLib::VecNormType norm_type) override
double getRelativeChangeFromPreviousTimestep(GlobalVector const &x, MathLib::VecNormType norm_type) override
double getRelativeChangeFromPreviousTimestep(GlobalVector const &x, MathLib::VecNormType norm_type) override
VecNormType
Norm type. Not declared as class type in order to use the members as integers.
Definition: LinAlgEnums.h:20
static NUMLIB_EXPORT VectorProvider & provider
void pushState(const double, GlobalVector const &x, InternalMatrixStorage const &) override
double getNewXWeight() const override
Returns .
double norm(MatrixOrVector const &x, MathLib::VecNormType type)
Definition: LinAlg.h:87
void getWeightedOldX(GlobalVector &y) const override
Returns .
std::unique_ptr< GlobalVector > _dx
Used to store .
#define OGS_FATAL(fmt,...)
Definition: Error.h:63
void axpy(MatrixOrVector &y, double const a, MatrixOrVector const &x)
Computes .
Definition: LinAlg.h:57
double computeRelativeChangeFromPreviousTimestep(GlobalVector const &x, GlobalVector const &x_old, MathLib::VecNormType norm_type)
void copy(MatrixOrVector const &x, MatrixOrVector &y)
Copies x to y.
Definition: LinAlg.h:36
double getRelativeChangeFromPreviousTimestep(GlobalVector const &x, MathLib::VecNormType norm_type) override
void scale(MatrixOrVector &x, double const a)
Scales x by a.
Definition: LinAlg.h:43
const double BDF_Coeffs[6][7]
Coefficients used in the backward differentiation formulas.