OGS 6.2.2-330-gf48c72f61.dirty.20200225212913
MatrixTranslator.cpp
Go to the documentation of this file.
1 
11 #include "MatrixTranslator.h"
12 
13 #include "MathLib/LinAlg/LinAlg.h"
14 
15 namespace NumLib
16 {
17 void MatrixTranslatorGeneral<ODESystemTag::FirstOrderImplicitQuasilinear>::
18  computeA(GlobalMatrix const& M, GlobalMatrix const& K,
19  GlobalMatrix& A) const
20 {
21  namespace LinAlg = MathLib::LinAlg;
22 
23  auto const dxdot_dx = _time_disc.getNewXWeight();
24 
25  // A = M * dxdot_dx + K
26  LinAlg::copy(M, A);
27  LinAlg::aypx(A, dxdot_dx, K);
28 }
29 
31  computeRhs(const GlobalMatrix& M, const GlobalMatrix& /*K*/,
32  const GlobalVector& b, GlobalVector& rhs) const
33 {
34  namespace LinAlg = MathLib::LinAlg;
35 
37  _time_disc.getWeightedOldX(tmp);
38 
39  // rhs = M * weighted_old_x + b
40  LinAlg::matMultAdd(M, tmp, b, rhs);
41 
43 }
44 
46  computeResidual(GlobalMatrix const& M, GlobalMatrix const& K,
47  GlobalVector const& b, GlobalVector const& x_new_timestep,
48  GlobalVector const& xdot, GlobalVector& res) const
49 {
50  namespace LinAlg = MathLib::LinAlg;
51 
52  auto const& x_curr = _time_disc.getCurrentX(x_new_timestep);
53 
54  // res = M * x_dot + K * x_curr - b
55  LinAlg::matMult(M, xdot, res); // the local vector x_dot seems to be
56  // necessary because of this multiplication
57  LinAlg::matMultAdd(K, x_curr, res, res);
58  LinAlg::axpy(res, -1.0, b);
59 }
60 
62  computeJacobian(GlobalMatrix const& Jac_in, GlobalMatrix& Jac_out) const
63 {
64  namespace LinAlg = MathLib::LinAlg;
65 
66  LinAlg::copy(Jac_in, Jac_out);
67 }
68 
70  computeA(GlobalMatrix const& M, GlobalMatrix const& /*K*/,
71  GlobalMatrix& A) const
72 {
73  namespace LinAlg = MathLib::LinAlg;
74 
75  auto const dxdot_dx = _fwd_euler.getNewXWeight();
76 
77  // A = M * dxdot_dx
78  LinAlg::copy(M, A);
79  LinAlg::scale(A, dxdot_dx);
80 }
81 
83  computeRhs(const GlobalMatrix& M, const GlobalMatrix& K,
84  const GlobalVector& b, GlobalVector& rhs) const
85 {
86  namespace LinAlg = MathLib::LinAlg;
87 
89  _fwd_euler.getWeightedOldX(tmp);
90 
91  auto const& x_old = _fwd_euler.getXOld();
92 
93  // rhs = b + M * weighted_old_x - K * x_old
94  LinAlg::matMult(K, x_old, rhs); // rhs = K * x_old
95  LinAlg::aypx(rhs, -1.0, b); // rhs = b - K * x_old
96  LinAlg::matMultAdd(M, tmp, rhs, rhs); // rhs += M * weighted_old_x
97 
99 }
100 
102  computeResidual(GlobalMatrix const& M, GlobalMatrix const& K,
103  GlobalVector const& b, GlobalVector const& x_new_timestep,
104  GlobalVector const& xdot, GlobalVector& res) const
105 {
106  namespace LinAlg = MathLib::LinAlg;
107 
108  auto const& x_curr = _fwd_euler.getCurrentX(x_new_timestep);
109 
110  // res = M * x_dot + K * x_curr - b
111  LinAlg::matMult(M, xdot, res);
112  LinAlg::matMultAdd(K, x_curr, res, res);
113  LinAlg::axpy(res, -1.0, b);
114 }
115 
117  computeJacobian(GlobalMatrix const& Jac_in, GlobalMatrix& Jac_out) const
118 {
119  namespace LinAlg = MathLib::LinAlg;
120 
121  LinAlg::copy(Jac_in, Jac_out);
122 }
123 
126  computeA(GlobalMatrix const& M, GlobalMatrix const& K,
127  GlobalMatrix& A) const
128 {
129  namespace LinAlg = MathLib::LinAlg;
130 
131  auto const dxdot_dx = _crank_nicolson.getNewXWeight();
132  auto const theta = _crank_nicolson.getTheta();
133 
134  // A = theta * (M * dxdot_dx + K) + dxdot_dx * _M_bar
135  LinAlg::copy(M, A);
136  LinAlg::aypx(A, dxdot_dx, K); // A = M * dxdot_dx + K
137 
138  LinAlg::scale(A, theta); // A *= theta
139  LinAlg::axpy(A, dxdot_dx, _M_bar); // A += dxdot_dx * _M_bar
140 }
141 
144  computeRhs(const GlobalMatrix& M, const GlobalMatrix& /*K*/,
145  const GlobalVector& b, GlobalVector& rhs) const
146 {
147  namespace LinAlg = MathLib::LinAlg;
148 
150  _crank_nicolson.getWeightedOldX(tmp);
151 
152  auto const theta = _crank_nicolson.getTheta();
153 
154  // rhs = theta * (b + M * weighted_old_x) + _M_bar * weighted_old_x -
155  // _b_bar;
156  LinAlg::matMultAdd(M, tmp, b, rhs); // rhs = b + M * weighted_old_x
157 
158  LinAlg::scale(rhs, theta); // rhs *= theta
159  LinAlg::matMultAdd(_M_bar, tmp, rhs, rhs); // rhs += _M_bar * weighted_old_x
160  LinAlg::axpy(rhs, -1.0, _b_bar); // rhs -= b
161 
163 }
164 
167  computeResidual(GlobalMatrix const& M, GlobalMatrix const& K,
168  GlobalVector const& b, GlobalVector const& x_new_timestep,
169  GlobalVector const& xdot, GlobalVector& res) const
170 {
171  namespace LinAlg = MathLib::LinAlg;
172 
173  auto const& x_curr = _crank_nicolson.getCurrentX(x_new_timestep);
174  auto const theta = _crank_nicolson.getTheta();
175 
176  // res = theta * (M * x_dot + K*x_curr - b) + _M_bar * x_dot + _b_bar
177  LinAlg::matMult(M, xdot, res); // res = M * x_dot
178  LinAlg::matMultAdd(K, x_curr, res, res); // res += K * x_curr
179  LinAlg::axpy(res, -1.0, b); // res = M * x_dot + K * x_curr - b
180 
181  LinAlg::aypx(res, theta, _b_bar); // res = res * theta + _b_bar
182  LinAlg::matMultAdd(_M_bar, xdot, res, res); // rs += _M_bar * x_dot
183 }
184 
187  computeJacobian(GlobalMatrix const& Jac_in, GlobalMatrix& Jac_out) const
188 {
189  namespace LinAlg = MathLib::LinAlg;
190 
191  auto const dxdot_dx = _crank_nicolson.getNewXWeight();
192  auto const theta = _crank_nicolson.getTheta();
193 
194  // J = theta * Jac + dxdot_dx * _M_bar
195  LinAlg::copy(Jac_in, Jac_out);
196  LinAlg::scale(Jac_out, theta);
197  LinAlg::axpy(Jac_out, dxdot_dx, _M_bar);
198 }
199 
202  pushMatrices(GlobalMatrix const& M, GlobalMatrix const& K,
203  GlobalVector const& b)
204 {
205  namespace LinAlg = MathLib::LinAlg;
206 
207  auto const theta = _crank_nicolson.getTheta();
208 
209  // Note: using x_old here is correct, since this method is called from
210  // within
211  // CrankNicolson::pushState() __after__ x_old has been updated to the
212  // result
213  // from the timestep just finished.
214  auto const& x_old = _crank_nicolson.getXOld();
215 
216  // _M_bar = (1.0-theta) * M;
217  LinAlg::copy(M, _M_bar);
218  LinAlg::scale(_M_bar, 1.0 - theta);
219 
220  // _b_bar = (1.0-theta) * (K * x_old - b)
221  LinAlg::matMult(K, x_old, _b_bar);
222  LinAlg::axpy(_b_bar, -1.0, b);
223  LinAlg::scale(_b_bar, 1.0 - theta);
224 }
225 
226 } // namespace NumLib
void aypx(MatrixOrVector &y, double const a, MatrixOrVector const &x)
Computes .
Definition: LinAlg.h:51
void matMult(Matrix const &A, Vector const &x, Vector &y)
Definition: LinAlg.h:114
virtual GlobalVector & getVector()=0
Get an uninitialized vector.
static NUMLIB_EXPORT VectorProvider & provider
virtual void releaseVector(GlobalVector const &x)=0
void axpy(MatrixOrVector &y, double const a, MatrixOrVector const &x)
Computes .
Definition: LinAlg.h:58
void copy(MatrixOrVector const &x, MatrixOrVector &y)
Copies x to y.
Definition: LinAlg.h:37
void matMultAdd(Matrix const &A, Vector const &v1, Vector const &v2, Vector &v3)
Definition: LinAlg.h:127
void scale(MatrixOrVector &x, double const a)
Scales x by a.
Definition: LinAlg.h:44