OGS 6.2.1-97-g73d1aeda3
CentralDifferencesJacobianAssembler.cpp
Go to the documentation of this file.
1 
11 #include "BaseLib/Error.h"
14 
15 namespace ProcessLib
16 {
18  std::vector<double>&& absolute_epsilons)
19  : _absolute_epsilons(std::move(absolute_epsilons))
20 {
21  if (_absolute_epsilons.empty())
22  {
23  OGS_FATAL("No values for the absolute epsilons have been given.");
24  }
25 }
26 
28  LocalAssemblerInterface& local_assembler, const double t,
29  const std::vector<double>& local_x_data,
30  const std::vector<double>& local_xdot_data, const double dxdot_dx,
31  const double dx_dx, std::vector<double>& local_M_data,
32  std::vector<double>& local_K_data, std::vector<double>& local_b_data,
33  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  OGS_FATAL(
38  "The number of specified epsilons (%u) and the number of local "
39  "d.o.f.s (%u) do not match, i.e., the latter is not divisable by "
40  "the former.",
41  _absolute_epsilons.size(), local_x_data.size());
42  }
43 
44  auto const num_r_c =
45  static_cast<Eigen::MatrixXd::Index>(local_x_data.size());
46 
47  auto const local_x =
48  MathLib::toVector<Eigen::VectorXd>(local_x_data, num_r_c);
49  auto const local_xdot =
50  MathLib::toVector<Eigen::VectorXd>(local_xdot_data, num_r_c);
51 
52  auto local_Jac = MathLib::createZeroedMatrix(local_Jac_data,
53  num_r_c, num_r_c);
54  _local_x_perturbed_data = local_x_data;
55 
56  auto const num_dofs_per_component =
57  local_x_data.size() / _absolute_epsilons.size();
58 
59  // Residual res := M xdot + K x - b
60  // Computing Jac := dres/dx
61  // = M dxdot/dx + dM/dx xdot + K dx/dx + dK/dx x - db/dx
62  // (Note: dM/dx and dK/dx actually have the second and
63  // third index transposed.)
64  // The loop computes the dM/dx, dK/dx and db/dx terms, the rest is computed
65  // afterwards.
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  _local_x_perturbed_data[i] += eps;
73  local_assembler.assemble(t, _local_x_perturbed_data, local_M_data,
74  local_K_data, local_b_data);
75 
76  _local_x_perturbed_data[i] = local_x_data[i] - eps;
79 
80  _local_x_perturbed_data[i] = local_x_data[i];
81 
82  if (!local_M_data.empty()) {
83  auto const local_M_p =
84  MathLib::toMatrix(local_M_data, num_r_c, num_r_c);
85  auto const local_M_m =
86  MathLib::toMatrix(_local_M_data, num_r_c, num_r_c);
87  local_Jac.col(i).noalias() +=
88  // dM/dxi * x_dot
89  (local_M_p - local_M_m) * local_xdot / (2.0 * eps);
90  local_M_data.clear();
91  _local_M_data.clear();
92  }
93  if (!local_K_data.empty()) {
94  auto const local_K_p =
95  MathLib::toMatrix(local_K_data, num_r_c, num_r_c);
96  auto const local_K_m =
97  MathLib::toMatrix(_local_K_data, num_r_c, num_r_c);
98  local_Jac.col(i).noalias() +=
99  // dK/dxi * x
100  (local_K_p - local_K_m) * local_x / (2.0 * eps);
101  local_K_data.clear();
102  _local_K_data.clear();
103  }
104  if (!local_b_data.empty()) {
105  auto const local_b_p =
106  MathLib::toVector<Eigen::VectorXd>(local_b_data, num_r_c);
107  auto const local_b_m =
108  MathLib::toVector<Eigen::VectorXd>(_local_b_data, num_r_c);
109  local_Jac.col(i).noalias() -=
110  // db/dxi
111  (local_b_p - local_b_m) / (2.0 * eps);
112  local_b_data.clear();
113  _local_b_data.clear();
114  }
115  }
116 
117  // Assemble with unperturbed local x.
118  local_assembler.assemble(t, local_x_data, local_M_data, local_K_data,
119  local_b_data);
120 
121  // Compute remaining terms of the Jacobian.
122  if (dxdot_dx != 0.0 && !local_M_data.empty()) {
123  auto local_M = MathLib::toMatrix(local_M_data, num_r_c, num_r_c);
124  local_Jac.noalias() += local_M * dxdot_dx;
125  }
126  if (dx_dx != 0.0 && !local_K_data.empty()) {
127  auto local_K = MathLib::toMatrix(local_K_data, num_r_c, num_r_c);
128  local_Jac.noalias() += local_K * dx_dx;
129  }
130 }
131 
132 std::unique_ptr<CentralDifferencesJacobianAssembler>
134 {
136  config.checkConfigParameter("type", "CentralDifferences");
137 
138  // TODO make non-optional.
140  auto rel_eps = config.getConfigParameterOptional<std::vector<double>>(
141  "relative_epsilons");
143  auto comp_mag = config.getConfigParameterOptional<std::vector<double>>(
144  "component_magnitudes");
145 
146  if (!!rel_eps != !!comp_mag) {
147  OGS_FATAL(
148  "Either both or none of <relative_epsilons> and "
149  "<component_magnitudes> have to be specified.");
150  }
151 
152  std::vector<double> abs_eps;
153 
154  if (rel_eps) {
155  if (rel_eps->size() != comp_mag->size()) {
156  OGS_FATAL(
157  "The numbers of components of <relative_epsilons> and "
158  "<component_magnitudes> do not match.");
159  }
160 
161  abs_eps.resize(rel_eps->size());
162  for (std::size_t i=0; i<rel_eps->size(); ++i) {
163  abs_eps[i] = (*rel_eps)[i] * (*comp_mag)[i];
164  }
165  } else {
166  // By default 1e-8 is used as epsilon for all components.
167  // TODO: remove this default value.
168  abs_eps.emplace_back(1e-8);
169  }
170 
171  return std::make_unique<CentralDifferencesJacobianAssembler>(
172  std::move(abs_eps));
173 }
174 
175 } // namespace ProcessLib
std::unique_ptr< CentralDifferencesJacobianAssembler > createCentralDifferencesJacobianAssembler(BaseLib::ConfigTree const &config)
Eigen::Map< Matrix > createZeroedMatrix(std::vector< double > &data, Eigen::MatrixXd::Index rows, Eigen::MatrixXd::Index cols)
Definition: EigenMapTools.h:31
void checkConfigParameter(std::string const &param, T const &value) const
Eigen::Map< const Matrix > toMatrix(std::vector< double > const &data, Eigen::MatrixXd::Index rows, Eigen::MatrixXd::Index cols)
Definition: EigenMapTools.h:71
boost::optional< T > getConfigParameterOptional(std::string const &param) const
#define OGS_FATAL(fmt,...)
Definition: Error.h:63
virtual void assemble(double const t, std::vector< double > const &local_x, std::vector< double > &local_M_data, std::vector< double > &local_K_data, std::vector< double > &local_b_data)
void assembleWithJacobian(LocalAssemblerInterface &local_assembler, double const t, std::vector< double > const &local_x, std::vector< double > const &local_xdot, const double dxdot_dx, const double dx_dx, 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
CentralDifferencesJacobianAssembler(std::vector< double > &&absolute_epsilons)