Loading web-font TeX/Math/Italic
OGS
ProcessLib::ForwardDifferencesJacobianAssembler Class Referencefinal

Detailed Description

Assembles the Jacobian matrix using forward differences.

Definition at line 17 of file ForwardDifferencesJacobianAssembler.h.

#include <ForwardDifferencesJacobianAssembler.h>

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

Public Member Functions

 ForwardDifferencesJacobianAssembler (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_b_data, std::vector< double > &local_Jac_data) override
 
std::unique_ptr< AbstractJacobianAssemblercopy () 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 > &)
 
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

◆ ForwardDifferencesJacobianAssembler()

ProcessLib::ForwardDifferencesJacobianAssembler::ForwardDifferencesJacobianAssembler ( 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 18 of file ForwardDifferencesJacobianAssembler.cpp.

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

References _absolute_epsilons, and OGS_FATAL.

Member Function Documentation

◆ assembleWithJacobian()

void ProcessLib::ForwardDifferencesJacobianAssembler::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_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.

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

Implements ProcessLib::AbstractJacobianAssembler.

Definition at line 28 of file ForwardDifferencesJacobianAssembler.cpp.

33{
34 std::vector<double> local_M_data(local_Jac_data.size());
35 std::vector<double> local_K_data(local_Jac_data.size());
36
37 // TODO do not check in every call.
38 if (local_x_data.size() % _absolute_epsilons.size() != 0)
39 {
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 x = MathLib::toVector<Eigen::VectorXd>(local_x_data, num_r_c);
51 auto const x_prev =
52 MathLib::toVector<Eigen::VectorXd>(local_x_prev_data, num_r_c);
53
54 Eigen::VectorXd const local_xdot = (x - x_prev) / dt;
55
56 auto local_Jac =
57 MathLib::createZeroedMatrix(local_Jac_data, num_r_c, num_r_c);
58
59 auto const num_dofs_per_component =
60 local_x_data.size() / _absolute_epsilons.size();
61
62 // Assemble with unperturbed local x to get M0, K0, and b0 used in the
63 // finite differences below.
64 local_assembler.assemble(t, dt, local_x_data, local_x_prev_data,
65 local_M_data, local_K_data, local_b_data);
66
67 auto const vds = local_assembler.getVectorDeformationSegment();
68
69 // Residual res := M xdot + K x - b
70 // Computing Jac := dres/dx
71 // = M dxdot/dx + dM/dx xdot + K dx/dx + dK/dx x - db/dx
72 // with dxdot/dx = 1/dt and dx/dx = 1
73 // (Note: dM/dx and dK/dx actually have the second and
74 // third index transposed.)
75 // The loop computes the dM/dx, dK/dx and db/dx terms, the rest is computed
76 // afterwards.
77 for (Eigen::MatrixXd::Index i = 0; i < num_r_c; ++i)
78 {
79 if ((vds != std::nullopt) && i >= vds->start_index &&
80 (i < (vds->start_index + vds->size)))
81 {
82 // Avoid to compute the analytic block
83 continue;
84 }
85
86 // assume that local_x_data is ordered by component.
87 auto const component = i / num_dofs_per_component;
88 auto const eps = _absolute_epsilons[component];
89
90 // Assemble with perturbed local x.
91 _local_x_perturbed_data = local_x_data;
92 _local_x_perturbed_data[i] = local_x_data[i] + eps;
93
94 local_assembler.assemble(t, dt, _local_x_perturbed_data,
95 local_x_prev_data, _local_M_data,
97
98 if (!local_M_data.empty() && !_local_M_data.empty())
99 {
100 auto const local_M_0 =
101 MathLib::toMatrix(local_M_data, num_r_c, num_r_c);
102 auto const local_M_p =
103 MathLib::toMatrix(_local_M_data, num_r_c, num_r_c);
104 local_Jac.col(i).noalias() +=
105 // dM/dxi * x_dot
106 (local_M_p - local_M_0) * local_xdot / eps;
107 _local_M_data.clear();
108 }
109 if (!local_K_data.empty() && !_local_K_data.empty())
110 {
111 auto const local_K_0 =
112 MathLib::toMatrix(local_K_data, num_r_c, num_r_c);
113 auto const local_K_p =
114 MathLib::toMatrix(_local_K_data, num_r_c, num_r_c);
115
116 local_Jac.col(i).noalias() +=
117 // dK/dxi * x
118 (local_K_p - local_K_0) * x / eps;
119 _local_K_data.clear();
120 }
121 if (!local_b_data.empty() && !_local_b_data.empty())
122 {
123 auto const local_b_0 =
124 MathLib::toVector<Eigen::VectorXd>(local_b_data, num_r_c);
125 auto const local_b_p =
127 local_Jac.col(i).noalias() -= (local_b_p - local_b_0) / eps;
128 _local_b_data.clear();
129 }
130 }
131
132 // Assemble with unperturbed local x, i.e. compute M dxdot/dx + K dx/dx =
133 // M/dt + K
134 // Compute remaining terms of the Jacobian.
135 if (!local_M_data.empty())
136 {
137 auto local_M = MathLib::toMatrix(local_M_data, num_r_c, num_r_c);
138 local_Jac.noalias() += local_M / dt;
139 }
140 if (!local_K_data.empty())
141 {
142 auto local_K = MathLib::toMatrix(local_K_data, num_r_c, num_r_c);
143 local_Jac.noalias() += local_K;
144 }
145
146 // Move the M and K contributions to the residuum for evaluation of nodal
147 // forces, flow rates, and the like. Cleaning up the M's and K's storage so
148 // it is not accounted for twice.
149 auto b = [&]()
150 {
151 if (!local_b_data.empty())
152 {
153 return MathLib::toVector<Eigen::VectorXd>(local_b_data, num_r_c);
154 }
156 num_r_c);
157 }();
158
159 if (!local_M_data.empty())
160 {
161 auto M = MathLib::toMatrix(local_M_data, num_r_c, num_r_c);
162 b -= M * (x - x_prev) / dt;
163 local_M_data.clear();
164 }
165 if (!local_K_data.empty())
166 {
167 auto K = MathLib::toMatrix(local_K_data, num_r_c, num_r_c);
168
169 // Note: The deformation segment of \c b is already computed as
170 // int{B^T sigma}dA, which is identical to K_uu * u. Therefore the
171 // corresponding K block is set to zero.
172 if (vds != std::nullopt)
173 {
174 K.block(vds->start_index, vds->start_index, vds->size, vds->size)
175 .setZero();
176 }
177
178 b -= K * x;
179 local_K_data.clear();
180 }
181}
Eigen::Map< Vector > createZeroedVector(std::vector< double > &data, Eigen::VectorXd::Index size)
Eigen::Map< const Vector > toVector(std::vector< double > const &data, Eigen::VectorXd::Index size)
Creates an Eigen mapped vector from the given data vector.
Eigen::Map< Matrix > createZeroedMatrix(std::vector< double > &data, Eigen::MatrixXd::Index rows, Eigen::MatrixXd::Index cols)
Eigen::Map< const Matrix > toMatrix(std::vector< double > const &data, Eigen::MatrixXd::Index rows, Eigen::MatrixXd::Index cols)

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

◆ copy()

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

Implements ProcessLib::AbstractJacobianAssembler.

Definition at line 184 of file ForwardDifferencesJacobianAssembler.cpp.

185{
186 return std::make_unique<ForwardDifferencesJacobianAssembler>(*this);
187}

Member Data Documentation

◆ _absolute_epsilons

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

◆ _local_b_data

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

Definition at line 62 of file ForwardDifferencesJacobianAssembler.h.

Referenced by assembleWithJacobian().

◆ _local_K_data

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

Definition at line 61 of file ForwardDifferencesJacobianAssembler.h.

Referenced by assembleWithJacobian().

◆ _local_M_data

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

Definition at line 60 of file ForwardDifferencesJacobianAssembler.h.

Referenced by assembleWithJacobian().

◆ _local_x_perturbed_data

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

Definition at line 63 of file ForwardDifferencesJacobianAssembler.h.

Referenced by assembleWithJacobian().


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