OGS
List of all members | Public Member Functions | Private Attributes
ProcessLib::CentralDifferencesJacobianAssembler Class Referencefinal

Detailed Description

Assembles the Jacobian matrix using central differences.

Definition at line 25 of file CentralDifferencesJacobianAssembler.h.

#include <CentralDifferencesJacobianAssembler.h>

Inheritance diagram for ProcessLib::CentralDifferencesJacobianAssembler:
[legend]
Collaboration diagram for ProcessLib::CentralDifferencesJacobianAssembler:
[legend]

Public Member Functions

 CentralDifferencesJacobianAssembler (std::vector< double > &&absolute_epsilons)
 
void assembleWithJacobian (LocalAssemblerInterface &local_assembler, double const t, double const dt, std::vector< double > const &local_x_data, std::vector< double > const &local_x_prev_data, std::vector< double > &local_M_data, std::vector< double > &local_K_data, std::vector< double > &local_b_data, std::vector< double > &local_Jac_data) override
 
std::unique_ptr< AbstractJacobianAssembler > copy () const override
 
- Public Member Functions inherited from ProcessLib::AbstractJacobianAssembler
virtual void assembleWithJacobianForStaggeredScheme (LocalAssemblerInterface &, double const, double const, Eigen::VectorXd const &, Eigen::VectorXd const &, int const, std::vector< double > &, std::vector< double > &, std::vector< double > &, std::vector< double > &)
 
virtual ~AbstractJacobianAssembler ()=default
 

Private Attributes

std::vector< double > const _absolute_epsilons
 
std::vector< double > _local_M_data
 
std::vector< double > _local_K_data
 
std::vector< double > _local_b_data
 
std::vector< double > _local_x_perturbed_data
 

Constructor & Destructor Documentation

◆ CentralDifferencesJacobianAssembler()

ProcessLib::CentralDifferencesJacobianAssembler::CentralDifferencesJacobianAssembler ( std::vector< double > && absolute_epsilons)
explicit

Constructs a new instance.

Parameters
absolute_epsilonsperturbations of the components of the local solution vector used for evaluating the finite differences.
Note
The size of absolute_epsilons defines the "number of components" of the local solution vector (This is not the number of elements of the vector!). Therefore the size of the local solution vector must be divisible by the size of absolute_epsilons. This is the only consistency check performed. It is not checked whether said "number of components" is sensible. E.g., one could pass one epsilon per node, which would be valid but would not make sense at all.

Definition at line 20 of file CentralDifferencesJacobianAssembler.cpp.

22 : _absolute_epsilons(std::move(absolute_epsilons))
23{
24 if (_absolute_epsilons.empty())
25 {
26 OGS_FATAL("No values for the absolute epsilons have been given.");
27 }
28}
OGS_FATAL
#define OGS_FATAL(...)
Definition Error.h:26

References _absolute_epsilons, and OGS_FATAL.

Member Function Documentation

◆ assembleWithJacobian()

void ProcessLib::CentralDifferencesJacobianAssembler::assembleWithJacobian ( LocalAssemblerInterface & local_assembler,
double const t,
double const dt,
std::vector< double > const & local_x_data,
std::vector< double > const & local_x_prev_data,
std::vector< double > & local_M_data,
std::vector< double > & local_K_data,
std::vector< double > & local_b_data,
std::vector< double > & local_Jac_data )
overridevirtual

Assembles the Jacobian, the matrices \(M\) and \(K\), and the vector \(b\). For the assembly the assemble() method of the given local_assembler is called several times and the Jacobian is built from finite differences. The number of calls of the assemble() method is \(2N+1\) if \(N\) is the size of local_x_data.

Attention
It is assumed that the local vectors and matrices are ordered by component.

Implements ProcessLib::AbstractJacobianAssembler.

Definition at line 30 of file CentralDifferencesJacobianAssembler.cpp.

36{
37 // TODO do not check in every call.
38 if (local_x_data.size() % _absolute_epsilons.size() != 0)
39 {
40 OGS_FATAL(
41 "The number of specified epsilons ({:d}) and the number of local "
42 "d.o.f.s ({:d}) do not match, i.e., the latter is not divisible by "
43 "the former.",
44 _absolute_epsilons.size(), local_x_data.size());
45 }
46
47 auto const num_r_c =
48 static_cast<Eigen::MatrixXd::Index>(local_x_data.size());
49
50 auto const local_x =
51 MathLib::toVector<Eigen::VectorXd>(local_x_data, num_r_c);
52 auto const local_x_prev =
53 MathLib::toVector<Eigen::VectorXd>(local_x_prev_data, num_r_c);
54 Eigen::VectorXd const local_xdot = (local_x - local_x_prev) / dt;
55
56 auto local_Jac =
57 MathLib::createZeroedMatrix(local_Jac_data, num_r_c, num_r_c);
58 _local_x_perturbed_data = local_x_data;
59
60 auto const num_dofs_per_component =
61 local_x_data.size() / _absolute_epsilons.size();
62
63 // Residual res := M xdot + K x - b
64 // Computing Jac := dres/dx
65 // = M dxdot/dx + dM/dx xdot + K dx/dx + dK/dx x - db/dx
66 // with dxdot/dx = 1/dt and dx/dx = 1
67 // (Note: dM/dx and dK/dx actually have the second and
68 // third index transposed.)
69 // The loop computes the dM/dx, dK/dx and db/dx terms, the rest is computed
70 // afterwards.
71 for (Eigen::MatrixXd::Index i = 0; i < num_r_c; ++i)
72 {
73 // assume that local_x_data is ordered by component.
74 auto const component = i / num_dofs_per_component;
75 auto const eps = _absolute_epsilons[component];
76
77 _local_x_perturbed_data[i] += eps;
78 local_assembler.assemble(t, dt, _local_x_perturbed_data,
79 local_x_prev_data, local_M_data, local_K_data,
80 local_b_data);
81
82 _local_x_perturbed_data[i] = local_x_data[i] - eps;
83 local_assembler.assemble(t, dt, _local_x_perturbed_data,
84 local_x_prev_data, _local_M_data,
85 _local_K_data, _local_b_data);
86
87 _local_x_perturbed_data[i] = local_x_data[i];
88
89 if (!local_M_data.empty())
90 {
91 auto const local_M_p =
92 MathLib::toMatrix(local_M_data, num_r_c, num_r_c);
93 auto const local_M_m =
94 MathLib::toMatrix(_local_M_data, num_r_c, num_r_c);
95 local_Jac.col(i).noalias() +=
96 // dM/dxi * x_dot
97 (local_M_p - local_M_m) * local_xdot / (2.0 * eps);
98 local_M_data.clear();
99 _local_M_data.clear();
100 }
101 if (!local_K_data.empty())
102 {
103 auto const local_K_p =
104 MathLib::toMatrix(local_K_data, num_r_c, num_r_c);
105 auto const local_K_m =
106 MathLib::toMatrix(_local_K_data, num_r_c, num_r_c);
107 local_Jac.col(i).noalias() +=
108 // dK/dxi * x
109 (local_K_p - local_K_m) * local_x / (2.0 * eps);
110 local_K_data.clear();
111 _local_K_data.clear();
112 }
113 if (!local_b_data.empty())
114 {
115 auto const local_b_p =
116 MathLib::toVector<Eigen::VectorXd>(local_b_data, num_r_c);
117 auto const local_b_m =
118 MathLib::toVector<Eigen::VectorXd>(_local_b_data, num_r_c);
119 local_Jac.col(i).noalias() -=
120 // db/dxi
121 (local_b_p - local_b_m) / (2.0 * eps);
122 local_b_data.clear();
123 _local_b_data.clear();
124 }
125 }
126
127 // Assemble with unperturbed local x.
128 local_assembler.assemble(t, dt, local_x_data, local_x_prev_data,
129 local_M_data, local_K_data, local_b_data);
130
131 // Compute remaining terms of the Jacobian.
132 if (!local_M_data.empty())
133 {
134 auto local_M = MathLib::toMatrix(local_M_data, num_r_c, num_r_c);
135 local_Jac.noalias() += local_M / dt;
136 }
137 if (!local_K_data.empty())
138 {
139 auto local_K = MathLib::toMatrix(local_K_data, num_r_c, num_r_c);
140 local_Jac.noalias() += local_K;
141 }
142
143 // Move the M and K contributions to the residuum for evaluation of nodal
144 // forces, flow rates, and the like. Cleaning up the M's and K's storage so
145 // it is not accounted for twice.
146 auto b = [&]()
147 {
148 if (!local_b_data.empty())
149 {
150 return MathLib::toVector<Eigen::VectorXd>(local_b_data, num_r_c);
151 }
152 return MathLib::createZeroedVector<Eigen::VectorXd>(local_b_data,
153 num_r_c);
154 }();
155
156 if (!local_M_data.empty())
157 {
158 auto M = MathLib::toMatrix(local_M_data, num_r_c, num_r_c);
159 b -= M * local_xdot;
160 local_M_data.clear();
161 }
162 if (!local_K_data.empty())
163 {
164 auto K = MathLib::toMatrix(local_K_data, num_r_c, num_r_c);
165 b -= K * local_x;
166 local_K_data.clear();
167 }
168}
MathLib::createZeroedMatrix
Eigen::Map< Matrix > createZeroedMatrix(std::vector< double > &data, Eigen::MatrixXd::Index rows, Eigen::MatrixXd::Index cols)
Definition EigenMapTools.h:32
MathLib::toMatrix
Eigen::Map< const Matrix > toMatrix(std::vector< double > const &data, Eigen::MatrixXd::Index rows, Eigen::MatrixXd::Index cols)
Definition EigenMapTools.h:72

References _absolute_epsilons, _local_b_data, _local_K_data, _local_M_data, _local_x_perturbed_data, ProcessLib::LocalAssemblerInterface::assemble(), MathLib::createZeroedMatrix(), OGS_FATAL, and MathLib::toMatrix().

◆ copy()

std::unique_ptr< AbstractJacobianAssembler > ProcessLib::CentralDifferencesJacobianAssembler::copy ( ) const
overridevirtual

Implements ProcessLib::AbstractJacobianAssembler.

Definition at line 220 of file CentralDifferencesJacobianAssembler.cpp.

221{
222 return std::make_unique<CentralDifferencesJacobianAssembler>(*this);
223}

Member Data Documentation

◆ _absolute_epsilons

std::vector<double> const ProcessLib::CentralDifferencesJacobianAssembler::_absolute_epsilons
private

Definition at line 66 of file CentralDifferencesJacobianAssembler.h.

Referenced by CentralDifferencesJacobianAssembler(), and assembleWithJacobian().

◆ _local_b_data

std::vector<double> ProcessLib::CentralDifferencesJacobianAssembler::_local_b_data
private

Definition at line 72 of file CentralDifferencesJacobianAssembler.h.

Referenced by assembleWithJacobian().

◆ _local_K_data

std::vector<double> ProcessLib::CentralDifferencesJacobianAssembler::_local_K_data
private

Definition at line 71 of file CentralDifferencesJacobianAssembler.h.

Referenced by assembleWithJacobian().

◆ _local_M_data

std::vector<double> ProcessLib::CentralDifferencesJacobianAssembler::_local_M_data
private

Definition at line 70 of file CentralDifferencesJacobianAssembler.h.

Referenced by assembleWithJacobian().

◆ _local_x_perturbed_data

std::vector<double> ProcessLib::CentralDifferencesJacobianAssembler::_local_x_perturbed_data
private

Definition at line 73 of file CentralDifferencesJacobianAssembler.h.

Referenced by assembleWithJacobian().

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