OGS
ForwardDifferencesJacobianAssembler.cpp
Go to the documentation of this file.
1
11
12#include "BaseLib/Error.h"
15
16namespace ProcessLib
17{
19 std::vector<double>&& absolute_epsilons)
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}
27
29 LocalAssemblerInterface& local_assembler, const double t, double const dt,
30 const std::vector<double>& local_x_data,
31 const std::vector<double>& local_x_prev_data,
32 std::vector<double>& local_M_data, std::vector<double>& local_K_data,
33 std::vector<double>& local_b_data, std::vector<double>& local_Jac_data)
34{
35 // TODO do not check in every call.
36 if (local_x_data.size() % _absolute_epsilons.size() != 0)
37 {
39 "The number of specified epsilons ({:d}) and the number of local "
40 "d.o.f.s ({:d}) do not match, i.e., the latter is not divisible by "
41 "the former.",
42 _absolute_epsilons.size(), local_x_data.size());
43 }
44
45 auto const num_r_c =
46 static_cast<Eigen::MatrixXd::Index>(local_x_data.size());
47
48 auto const x = MathLib::toVector<Eigen::VectorXd>(local_x_data, num_r_c);
49 auto const x_prev =
50 MathLib::toVector<Eigen::VectorXd>(local_x_prev_data, num_r_c);
51
52 auto local_Jac =
53 MathLib::createZeroedMatrix(local_Jac_data, num_r_c, num_r_c);
54
55 auto const num_dofs_per_component =
56 local_x_data.size() / _absolute_epsilons.size();
57
58 // Assemble with unperturbed local x to get M0, K0, and b0 used in the
59 // finite differences below.
60 local_assembler.assemble(t, dt, local_x_data, local_x_prev_data,
61 local_M_data, local_K_data, local_b_data);
62
63 // Residual res := M xdot + K x - b
64 // Computing Jac := dres/dx
65 // = d(M xdot)/dx + d(K x)/dx - db/dx
66 for (Eigen::MatrixXd::Index i = 0; i < num_r_c; ++i)
67 {
68 // assume that local_x_data is ordered by component.
69 auto const component = i / num_dofs_per_component;
70 auto const eps = _absolute_epsilons[component];
71
72 // Assemble with perturbed local x.
73 _local_x_perturbed_data = local_x_data;
74 _local_x_perturbed_data[i] = local_x_data[i] + eps;
75 auto const x_p = MathLib::toVector<Eigen::VectorXd>(
77
78 local_assembler.assemble(t, dt, _local_x_perturbed_data,
79 local_x_prev_data, _local_M_data,
81
82 if (!local_M_data.empty() && !_local_M_data.empty())
83 {
84 auto const local_M_0 =
85 MathLib::toMatrix(local_M_data, num_r_c, num_r_c);
86 auto const local_M_p =
87 MathLib::toMatrix(_local_M_data, num_r_c, num_r_c);
88 local_Jac.col(i).noalias() +=
89 (local_M_p * (x_p - x_prev) - local_M_0 * (x - x_prev)) /
90 (dt * eps);
91 _local_M_data.clear();
92 }
93 if (!local_K_data.empty() && !_local_K_data.empty())
94 {
95 auto const local_K_0 =
96 MathLib::toMatrix(local_K_data, num_r_c, num_r_c);
97 auto const local_K_p =
98 MathLib::toMatrix(_local_K_data, num_r_c, num_r_c);
99 local_Jac.col(i).noalias() +=
100 (local_K_p * x_p - local_K_0 * x) / eps;
101 _local_K_data.clear();
102 }
103 if (!local_b_data.empty() && !_local_b_data.empty())
104 {
105 auto const local_b_0 =
106 MathLib::toVector<Eigen::VectorXd>(local_b_data, num_r_c);
107 auto const local_b_p =
108 MathLib::toVector<Eigen::VectorXd>(_local_b_data, num_r_c);
109 local_Jac.col(i).noalias() -= (local_b_p - local_b_0) / eps;
110 _local_b_data.clear();
111 }
112 }
113
114 local_M_data.clear();
115 local_K_data.clear();
116 local_b_data.clear();
117 // Assemble with unperturbed local x to keep the internal assembler state
118 // for next iteration or time step.
119 local_assembler.assemble(t, dt, local_x_data, local_x_prev_data,
120 local_M_data, local_K_data, local_b_data);
121
122 // Move the M and K contributions to the residuum for evaluation of nodal
123 // forces, flow rates, and the like. Cleaning up the M's and K's storage so
124 // it is not accounted for twice.
125 auto b = [&]()
126 {
127 if (!local_b_data.empty())
128 {
129 return MathLib::toVector<Eigen::VectorXd>(local_b_data, num_r_c);
130 }
131 return MathLib::createZeroedVector<Eigen::VectorXd>(local_b_data,
132 num_r_c);
133 }();
134
135 if (!local_M_data.empty())
136 {
137 auto M = MathLib::toMatrix(local_M_data, num_r_c, num_r_c);
138 b -= M * (x - x_prev) / dt;
139 local_M_data.clear();
140 }
141 if (!local_K_data.empty())
142 {
143 auto K = MathLib::toMatrix(local_K_data, num_r_c, num_r_c);
144 b -= K * x;
145 local_K_data.clear();
146 }
147}
148
149std::unique_ptr<AbstractJacobianAssembler>
151{
152 return std::make_unique<ForwardDifferencesJacobianAssembler>(*this);
153}
154
155} // namespace ProcessLib
#define OGS_FATAL(...)
Definition Error.h:26
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_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
virtual void assemble(double const t, double const dt, std::vector< double > const &local_x, std::vector< double > const &local_x_prev, std::vector< double > &local_M_data, std::vector< double > &local_K_data, std::vector< double > &local_b_data)
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)