OGS 6.2.1-97-g73d1aeda3
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  const double /*t*/,
83  GlobalVector const& x,
84  InternalMatrixStorage const& /*strg*/)
85 {
86  namespace LinAlg = MathLib::LinAlg;
87  // TODO use boost circular buffer?
88 
89  // until _xs_old is filled, lower-order BDF formulas are used.
90  if (_xs_old.size() < _num_steps)
91  {
92  _xs_old.push_back(&NumLib::GlobalVectorProvider::provider.getVector(x));
93  }
94  else
95  {
96  LinAlg::copy(x, *_xs_old[_offset]);
97  _offset =
98  (_offset + 1) % _num_steps; // treat _xs_old as a circular buffer
99  }
100 }
101 
102 namespace detail
103 {
105 const double BDF_Coeffs[6][7] = {
106  // leftmost column: weight of the solution at the new timestep
107  // signs of columns > 1 are flipped compared to standard BDF tableaus
108  {1.0, 1.0},
109  {1.5, 2.0, -0.5},
110  {11.0 / 6.0, 3.0, -1.5, 1.0 / 3.0},
111  {25.0 / 12.0, 4.0, -3.0, 4.0 / 3.0, -0.25},
112  {137.0 / 60.0, 5.0, -5.0, 10.0 / 3.0, -1.25, 1.0 / 5.0},
113  {147.0 / 60.0, 6.0, -7.5, 20.0 / 3.0, -3.75, 6.0 / 5.0, -1.0 / 6.0}
114  // coefficient of (for BDF(6), the oldest state, x_n, is always rightmost)
115  // x_+6, x_+5, x_+4, x_+3, x_+2, x_+1, x_n
116 };
117 } // namespace detail
118 
120 {
121  auto const k = eff_num_steps();
122  return detail::BDF_Coeffs[k - 1][0] / _delta_t;
123 }
124 
126 {
127  namespace LinAlg = MathLib::LinAlg;
128 
129  auto const k = eff_num_steps();
130  auto const* const BDFk = detail::BDF_Coeffs[k - 1];
131 
132  // compute linear combination \sum_{i=0}^{k-1} BDFk_{k-i} \cdot x_{n+i}
133  LinAlg::copy(*_xs_old[_offset], y); // _xs_old[offset] = x_n
134  LinAlg::scale(y, BDFk[k]);
135 
136  for (unsigned i = 1; i < k; ++i)
137  {
138  auto const off = (_offset + i) % k;
139  LinAlg::axpy(y, BDFk[k - i], *_xs_old[off]);
140  }
141 
142  LinAlg::scale(y, 1.0 / _delta_t);
143 }
144 
145 } // 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.