OGS 6.2.0-97-g4a610c866
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 210 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 125 of file MatrixTranslator.cpp.

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

127 {
128  namespace LinAlg = MathLib::LinAlg;
129 
130  auto const dxdot_dx = _crank_nicolson.getNewXWeight();
131  auto const theta = _crank_nicolson.getTheta();
132 
133  // A = theta * (M * dxdot_dx + K) + dxdot_dx * _M_bar
134  LinAlg::copy(M, A);
135  LinAlg::aypx(A, dxdot_dx, K); // A = M * dxdot_dx + K
136 
137  LinAlg::scale(A, theta); // A *= theta
138  LinAlg::axpy(A, dxdot_dx, _M_bar); // A += dxdot_dx * _M_bar
139 }
void aypx(MatrixOrVector &y, double const a, MatrixOrVector const &x)
Computes .
Definition: LinAlg.h:50
double getNewXWeight() const override
Returns .
double getTheta() const
Returns .
void axpy(MatrixOrVector &y, double const a, MatrixOrVector const &x)
Computes .
Definition: LinAlg.h:57
void copy(MatrixOrVector const &x, MatrixOrVector &y)
Copies x to y.
Definition: LinAlg.h:36
void scale(MatrixOrVector &x, double const a)
Scales x by a.
Definition: LinAlg.h:43

◆ 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 186 of file MatrixTranslator.cpp.

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

187 {
188  namespace LinAlg = MathLib::LinAlg;
189 
190  auto const dxdot_dx = _crank_nicolson.getNewXWeight();
191  auto const theta = _crank_nicolson.getTheta();
192 
193  // J = theta * Jac + dxdot_dx * _M_bar
194  LinAlg::copy(Jac_in, Jac_out);
195  LinAlg::scale(Jac_out, theta);
196  LinAlg::axpy(Jac_out, dxdot_dx, _M_bar);
197 }
double getNewXWeight() const override
Returns .
double getTheta() const
Returns .
void axpy(MatrixOrVector &y, double const a, MatrixOrVector const &x)
Computes .
Definition: LinAlg.h:57
void copy(MatrixOrVector const &x, MatrixOrVector &y)
Copies x to y.
Definition: LinAlg.h:36
void scale(MatrixOrVector &x, double const a)
Scales x by a.
Definition: LinAlg.h:43

◆ 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 166 of file MatrixTranslator.cpp.

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

169 {
170  namespace LinAlg = MathLib::LinAlg;
171 
172  auto const& x_curr = _crank_nicolson.getCurrentX(x_new_timestep);
173  auto const theta = _crank_nicolson.getTheta();
174 
175  // res = theta * (M * x_dot + K*x_curr - b) + _M_bar * x_dot + _b_bar
176  LinAlg::matMult(M, xdot, res); // res = M * x_dot
177  LinAlg::matMultAdd(K, x_curr, res, res); // res += K * x_curr
178  LinAlg::axpy(res, -1.0, b); // res = M * x_dot + K * x_curr - b
179 
180  LinAlg::aypx(res, theta, _b_bar); // res = res * theta + _b_bar
181  LinAlg::matMultAdd(_M_bar, xdot, res, res); // rs += _M_bar * x_dot
182 }
void aypx(MatrixOrVector &y, double const a, MatrixOrVector const &x)
Computes .
Definition: LinAlg.h:50
void matMult(Matrix const &A, Vector const &x, Vector &y)
Definition: LinAlg.h:113
double getTheta() const
Returns .
void axpy(MatrixOrVector &y, double const a, MatrixOrVector const &x)
Computes .
Definition: LinAlg.h:57
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:126

◆ 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 143 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().

145 {
146  namespace LinAlg = MathLib::LinAlg;
147 
150 
151  auto const theta = _crank_nicolson.getTheta();
152 
153  // rhs = theta * (b + M * weighted_old_x) + _M_bar * weighted_old_x -
154  // _b_bar;
155  LinAlg::matMultAdd(M, tmp, b, rhs); // rhs = b + M * weighted_old_x
156 
157  LinAlg::scale(rhs, theta); // rhs *= theta
158  LinAlg::matMultAdd(_M_bar, tmp, rhs, rhs); // rhs += _M_bar * weighted_old_x
159  LinAlg::axpy(rhs, -1.0, _b_bar); // rhs -= b
160 
162 }
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:57
void matMultAdd(Matrix const &A, Vector const &v1, Vector const &v2, Vector &v3)
Definition: LinAlg.h:126
void scale(MatrixOrVector &x, double const a)
Scales x by a.
Definition: LinAlg.h:43

◆ 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 201 of file MatrixTranslator.cpp.

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

203 {
204  namespace LinAlg = MathLib::LinAlg;
205 
206  auto const theta = _crank_nicolson.getTheta();
207 
208  // Note: using x_old here is correct, since this method is called from
209  // within
210  // CrankNicolson::pushState() __after__ x_old has been updated to the
211  // result
212  // from the timestep just finished.
213  auto const& x_old = _crank_nicolson.getXOld();
214 
215  // _M_bar = (1.0-theta) * M;
216  LinAlg::copy(M, _M_bar);
217  LinAlg::scale(_M_bar, 1.0 - theta);
218 
219  // _b_bar = (1.0-theta) * (K * x_old - b)
220  LinAlg::matMult(K, x_old, _b_bar);
221  LinAlg::axpy(_b_bar, -1.0, b);
222  LinAlg::scale(_b_bar, 1.0 - theta);
223 }
void matMult(Matrix const &A, Vector const &x, Vector &y)
Definition: LinAlg.h:113
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:57
void copy(MatrixOrVector const &x, MatrixOrVector &y)
Copies x to y.
Definition: LinAlg.h:36
void scale(MatrixOrVector &x, double const a)
Scales x by a.
Definition: LinAlg.h:43

Member Data Documentation

◆ _b_bar

Used to adjust vectors assembled by the ODE.

See also
pushMatrices()

Definition at line 280 of file MatrixTranslator.h.

◆ _crank_nicolson

◆ _M_bar

Used to adjust matrices and vectors assembled by the ODE.

See also
pushMatrices()

Definition at line 278 of file MatrixTranslator.h.

◆ _tmp_id

ID of the vector storing intermediate computations.

Definition at line 284 of file MatrixTranslator.h.


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