OGS 6.1.0-1721-g6382411ad
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)
30 
31  auto& dx = *_dx;
32  MathLib::LinAlg::copy(x, dx); // copy x to dx.
33 
34  // dx = x - x_old --> x - dx --> dx
35  MathLib::LinAlg::axpy(dx, -1.0, x_old);
36  const double norm_dx = MathLib::LinAlg::norm(dx, norm_type);
37 
38  const double norm_x = MathLib::LinAlg::norm(x, norm_type);
39  if (norm_x > std::numeric_limits<double>::epsilon())
40  return norm_dx / norm_x;
41 
42  // Both of norm_x and norm_dx are close to zero
43  if (norm_dx < std::numeric_limits<double>::epsilon())
44  return 1.0;
45 
46  // Only norm_x is close to zero
47  return norm_dx / std::numeric_limits<double>::epsilon();
48 }
49 
51  GlobalVector const& x, MathLib::VecNormType norm_type)
52 {
53  return computeRelativeChangeFromPreviousTimestep(x, _x_old, norm_type);
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 {
72  x, *_xs_old[_offset], norm_type);
73 }
74 
76  GlobalVector const& x,
77  InternalMatrixStorage const&)
78 {
79  namespace LinAlg = MathLib::LinAlg;
80  // TODO use boost circular buffer?
81 
82  // until _xs_old is filled, lower-order BDF formulas are used.
83  if (_xs_old.size() < _num_steps)
84  {
85  _xs_old.push_back(&NumLib::GlobalVectorProvider::provider.getVector(x));
86  }
87  else
88  {
89  LinAlg::copy(x, *_xs_old[_offset]);
90  _offset =
91  (_offset + 1) % _num_steps; // treat _xs_old as a circular buffer
92  }
93 }
94 
95 namespace detail
96 {
98 const double BDF_Coeffs[6][7] = {
99  // leftmost column: weight of the solution at the new timestep
100  // signs of columns > 1 are flipped compared to standard BDF tableaus
101  {1.0, 1.0},
102  {1.5, 2.0, -0.5},
103  {11.0 / 6.0, 3.0, -1.5, 1.0 / 3.0},
104  {25.0 / 12.0, 4.0, -3.0, 4.0 / 3.0, -0.25},
105  {137.0 / 60.0, 5.0, -5.0, 10.0 / 3.0, -1.25, 1.0 / 5.0},
106  {147.0 / 60.0, 6.0, -7.5, 20.0 / 3.0, -3.75, 6.0 / 5.0, -1.0 / 6.0}
107  // coefficient of (for BDF(6), the oldest state, x_n, is always rightmost)
108  // x_+6, x_+5, x_+4, x_+3, x_+2, x_+1, x_n
109 };
110 }
111 
113 {
114  auto const k = eff_num_steps();
115  return detail::BDF_Coeffs[k - 1][0] / _delta_t;
116 }
117 
119 {
120  namespace LinAlg = MathLib::LinAlg;
121 
122  auto const k = eff_num_steps();
123  auto const* const BDFk = detail::BDF_Coeffs[k - 1];
124 
125  // compute linear combination \sum_{i=0}^{k-1} BDFk_{k-i} \cdot x_{n+i}
126  LinAlg::copy(*_xs_old[_offset], y); // _xs_old[offset] = x_n
127  LinAlg::scale(y, BDFk[k]);
128 
129  for (unsigned i = 1; i < k; ++i)
130  {
131  auto const off = (_offset + i) % k;
132  LinAlg::axpy(y, BDFk[k - i], *_xs_old[off]);
133  }
134 
135  LinAlg::scale(y, 1.0 / _delta_t);
136 }
137 
138 } // 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:71
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.