OGS
CentralDifferencesJacobianAssembler.cpp
Go to the documentation of this file.
1
12
13#include "BaseLib/ConfigTree.h"
14#include "BaseLib/Error.h"
17
18namespace ProcessLib
19{
21 std::vector<double>&& absolute_epsilons)
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}
29
31 LocalAssemblerInterface& local_assembler, const double t, double const dt,
32 const std::vector<double>& local_x_data,
33 const std::vector<double>& local_x_prev_data,
34 std::vector<double>& local_M_data, std::vector<double>& local_K_data,
35 std::vector<double>& local_b_data, std::vector<double>& local_Jac_data)
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 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
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,
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}
169
170std::unique_ptr<CentralDifferencesJacobianAssembler>
172{
174 config.checkConfigParameter("type", "CentralDifferences");
175
176 // TODO make non-optional.
178 auto rel_eps = config.getConfigParameterOptional<std::vector<double>>(
179 "relative_epsilons");
181 auto comp_mag = config.getConfigParameterOptional<std::vector<double>>(
182 "component_magnitudes");
183
184 if (!!rel_eps != !!comp_mag)
185 {
186 OGS_FATAL(
187 "Either both or none of <relative_epsilons> and "
188 "<component_magnitudes> have to be specified.");
189 }
190
191 std::vector<double> abs_eps;
192
193 if (rel_eps)
194 {
195 if (rel_eps->size() != comp_mag->size())
196 {
197 OGS_FATAL(
198 "The numbers of components of <relative_epsilons> and "
199 "<component_magnitudes> do not match.");
200 }
201
202 abs_eps.resize(rel_eps->size());
203 for (std::size_t i = 0; i < rel_eps->size(); ++i)
204 {
205 abs_eps[i] = (*rel_eps)[i] * (*comp_mag)[i];
206 }
207 }
208 else
209 {
210 // By default 1e-8 is used as epsilon for all components.
211 // TODO: remove this default value.
212 abs_eps.emplace_back(1e-8);
213 }
214
215 return std::make_unique<CentralDifferencesJacobianAssembler>(
216 std::move(abs_eps));
217}
218
219std::unique_ptr<AbstractJacobianAssembler>
221{
222 return std::make_unique<CentralDifferencesJacobianAssembler>(*this);
223}
224
225} // namespace ProcessLib
#define OGS_FATAL(...)
Definition Error.h:26
std::optional< T > getConfigParameterOptional(std::string const &param) const
void checkConfigParameter(std::string const &param, std::string_view const value) const
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
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)
std::unique_ptr< CentralDifferencesJacobianAssembler > createCentralDifferencesJacobianAssembler(BaseLib::ConfigTree const &config)