OGS 6.2.1-720-gf9530c21a.dirty.20191210150148
NumLib::MatrixTranslatorCrankNicolson< ODESystemTag::FirstOrderImplicitQuasilinear > Class Template Reference

Detailed Description

template<>
class NumLib::MatrixTranslatorCrankNicolson< ODESystemTag::FirstOrderImplicitQuasilinear >

GlobalMatrix translator for first order implicit quasi-linear ODEs, used with the CrankNicolson scheme.

See also
ODESystemTag::FirstOrderImplicitQuasilinear
Remarks
You might also want read the remarks on time discretization.

Definition at line 211 of file MatrixTranslator.h.

#include <MatrixTranslator.h>

Inheritance diagram for NumLib::MatrixTranslatorCrankNicolson< ODESystemTag::FirstOrderImplicitQuasilinear >:
Collaboration diagram for NumLib::MatrixTranslatorCrankNicolson< ODESystemTag::FirstOrderImplicitQuasilinear >:

Public Member Functions

 MatrixTranslatorCrankNicolson (CrankNicolson const &timeDisc)
 
 ~MatrixTranslatorCrankNicolson () override
 
void computeA (GlobalMatrix const &M, GlobalMatrix const &K, GlobalMatrix &A) const override
 
void computeRhs (const GlobalMatrix &M, const GlobalMatrix &, const GlobalVector &b, GlobalVector &rhs) const override
 
void computeResidual (GlobalMatrix const &M, GlobalMatrix const &K, GlobalVector const &b, GlobalVector const &x_new_timestep, GlobalVector const &xdot, GlobalVector &res) const override
 
void computeJacobian (GlobalMatrix const &Jac_in, GlobalMatrix &Jac_out) const override
 
void pushMatrices (GlobalMatrix const &M, GlobalMatrix const &K, GlobalVector const &b) override
 
- Public Member Functions inherited from NumLib::MatrixTranslator< ODESystemTag::FirstOrderImplicitQuasilinear >
virtual ~MatrixTranslator ()=default
 

Private Attributes

CrankNicolson const & _crank_nicolson
 
GlobalMatrix & _M_bar
 
GlobalVector & _b_bar
 
std::size_t _tmp_id = 0u
 ID of the vector storing intermediate computations. More...
 

Constructor & Destructor Documentation

◆ MatrixTranslatorCrankNicolson()

◆ ~MatrixTranslatorCrankNicolson()

Member Function Documentation

◆ computeA()

void NumLib::MatrixTranslatorCrankNicolson< ODESystemTag::FirstOrderImplicitQuasilinear >::computeA ( GlobalMatrix const &  M,
GlobalMatrix const &  K,
GlobalMatrix &  A 
) const
overridevirtual

Computes $ A = \theta \cdot (M \cdot \alpha + K) + \bar M \cdot \alpha //! $.

Implements NumLib::MatrixTranslator< ODESystemTag::FirstOrderImplicitQuasilinear >.

Definition at line 126 of file MatrixTranslator.cpp.

References MathLib::LinAlg::axpy(), MathLib::LinAlg::aypx(), MathLib::LinAlg::copy(), NumLib::FirstOrderImplicitQuasilinear, and MathLib::LinAlg::scale().

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 }
void aypx(MatrixOrVector &y, double const a, MatrixOrVector const &x)
Computes .
Definition: LinAlg.h:51
double getNewXWeight() const override
Returns .
double getTheta() const
Returns .
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 scale(MatrixOrVector &x, double const a)
Scales x by a.
Definition: LinAlg.h:44

◆ computeJacobian()

void NumLib::MatrixTranslatorCrankNicolson< ODESystemTag::FirstOrderImplicitQuasilinear >::computeJacobian ( GlobalMatrix const &  Jac_in,
GlobalMatrix &  Jac_out 
) const
overridevirtual

Computes $ \mathtt{Jac\_out} = \theta \cdot \mathtt{Jac\_in} + \bar M \cdot \alpha $.

Where Jac_in is the Jacobian as assembled by the ODE system, i.e. in the same fashion as for the BackwardEuler scheme.

Implements NumLib::MatrixTranslator< ODESystemTag::FirstOrderImplicitQuasilinear >.

Definition at line 187 of file MatrixTranslator.cpp.

References MathLib::LinAlg::axpy(), MathLib::LinAlg::copy(), NumLib::FirstOrderImplicitQuasilinear, and MathLib::LinAlg::scale().

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 }
double getNewXWeight() const override
Returns .
double getTheta() const
Returns .
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 scale(MatrixOrVector &x, double const a)
Scales x by a.
Definition: LinAlg.h:44

◆ computeResidual()

void NumLib::MatrixTranslatorCrankNicolson< ODESystemTag::FirstOrderImplicitQuasilinear >::computeResidual ( GlobalMatrix const &  M,
GlobalMatrix const &  K,
GlobalVector const &  b,
GlobalVector const &  x_new_timestep,
GlobalVector const &  xdot,
GlobalVector &  res 
) const
overridevirtual

Computes $ r = \theta \cdot (M \cdot \hat x + K \cdot x_C - b) + \bar //! M \cdot \hat x + \bar b $.

Implements NumLib::MatrixTranslator< ODESystemTag::FirstOrderImplicitQuasilinear >.

Definition at line 167 of file MatrixTranslator.cpp.

References MathLib::LinAlg::axpy(), MathLib::LinAlg::aypx(), NumLib::FirstOrderImplicitQuasilinear, MathLib::LinAlg::matMult(), and MathLib::LinAlg::matMultAdd().

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 }
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
double getTheta() const
Returns .
void axpy(MatrixOrVector &y, double const a, MatrixOrVector const &x)
Computes .
Definition: LinAlg.h:58
virtual GlobalVector const & getCurrentX(GlobalVector const &x_at_new_timestep) const
void matMultAdd(Matrix const &A, Vector const &v1, Vector const &v2, Vector &v3)
Definition: LinAlg.h:127

◆ computeRhs()

void NumLib::MatrixTranslatorCrankNicolson< ODESystemTag::FirstOrderImplicitQuasilinear >::computeRhs ( const GlobalMatrix &  M,
const GlobalMatrix &  ,
const GlobalVector &  b,
GlobalVector &  rhs 
) const
overridevirtual

Computes $ \mathtt{rhs} = \theta \cdot (M \cdot x_O + b) + \bar M //! \cdot x_O - \bar b $.

Implements NumLib::MatrixTranslator< ODESystemTag::FirstOrderImplicitQuasilinear >.

Definition at line 144 of file MatrixTranslator.cpp.

References MathLib::LinAlg::axpy(), NumLib::FirstOrderImplicitQuasilinear, NumLib::VectorProvider::getVector(), MathLib::LinAlg::matMultAdd(), NumLib::GlobalVectorProvider::provider, NumLib::VectorProvider::releaseVector(), and MathLib::LinAlg::scale().

146 {
147  namespace LinAlg = MathLib::LinAlg;
148 
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 }
virtual GlobalVector & getVector()=0
Get an uninitialized vector.
void getWeightedOldX(GlobalVector &y) const override
Returns .
static NUMLIB_EXPORT VectorProvider & provider
double getTheta() const
Returns .
std::size_t _tmp_id
ID of the vector storing intermediate computations.
virtual void releaseVector(GlobalVector const &x)=0
void axpy(MatrixOrVector &y, double const a, MatrixOrVector const &x)
Computes .
Definition: LinAlg.h:58
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

◆ pushMatrices()

void NumLib::MatrixTranslatorCrankNicolson< ODESystemTag::FirstOrderImplicitQuasilinear >::pushMatrices ( GlobalMatrix const &  M,
GlobalMatrix const &  K,
GlobalVector const &  b 
)
overridevirtual

Saves internal state for use in the successive timestep; computes $ \bar M $ and $ \bar b $.

$ \bar M $ and $ \bar b $ are computed as follows:

\begin{align} \bar M &= (1-\theta) \cdot M \\ \bar b &= (1-\theta) \cdot ( K \cdot x_n - b ) \end{align}

Where $ x_n $ is the solution at the timestep just finished.

Reimplemented from NumLib::MatrixTranslator< ODESystemTag::FirstOrderImplicitQuasilinear >.

Definition at line 202 of file MatrixTranslator.cpp.

References MathLib::LinAlg::axpy(), MathLib::LinAlg::copy(), MathLib::LinAlg::matMult(), and MathLib::LinAlg::scale().

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 }
void matMult(Matrix const &A, Vector const &x, Vector &y)
Definition: LinAlg.h:114
GlobalVector const & getXOld() const
Returns the solution from the preceding timestep.
double getTheta() const
Returns .
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 scale(MatrixOrVector &x, double const a)
Scales x by a.
Definition: LinAlg.h:44

Member Data Documentation

◆ _b_bar

Used to adjust vectors assembled by the ODE.

See also
pushMatrices()

Definition at line 281 of file MatrixTranslator.h.

◆ _crank_nicolson

◆ _M_bar

Used to adjust matrices and vectors assembled by the ODE.

See also
pushMatrices()

Definition at line 279 of file MatrixTranslator.h.

◆ _tmp_id

ID of the vector storing intermediate computations.

Definition at line 285 of file MatrixTranslator.h.


The documentation for this class was generated from the following files: