OGS 6.2.0-97-g4a610c866
CompareJacobiansJacobianAssembler.cpp
Go to the documentation of this file.
1 
11 
12 #include <sstream>
13 
15 
17 
18 namespace
19 {
21 void dump_py(std::ostream& fh, std::string const& var, double const val)
22 {
23  fh << var << " = " << val << '\n';
24 }
25 
27 void dump_py(std::ostream& fh, std::string const& var, std::size_t const val)
28 {
29  fh << var << " = " << val << '\n';
30 }
31 
33 template <typename Vec>
34 void dump_py_vec(std::ostream& fh, std::string const& var, Vec const& val)
35 {
36  fh << var << " = np.array([";
37  for (decltype(val.size()) i = 0; i < val.size(); ++i)
38  {
39  if (i != 0)
40  {
41  if (i % 8 == 0)
42  {
43  // Print at most eight entries on one line,
44  // indent with four spaces.
45  fh << ",\n ";
46  }
47  else
48  {
49  fh << ", ";
50  }
51  }
52  fh << val[i];
53  }
54  fh << "])\n";
55 }
56 
58 void dump_py(std::ostream& fh, std::string const& var,
59  std::vector<double> const& val)
60 {
61  dump_py_vec(fh, var, val);
62 }
63 
65 template <typename Derived>
66 void dump_py(std::ostream& fh, std::string const& var,
67  Eigen::ArrayBase<Derived> const& val,
68  std::integral_constant<int, 1>)
69 {
70  dump_py_vec(fh, var, val);
71 }
72 
74 template <typename Derived, int ColsAtCompileTime>
75 void dump_py(std::ostream& fh, std::string const& var,
76  Eigen::ArrayBase<Derived> const& val,
77  std::integral_constant<int, ColsAtCompileTime>)
78 {
79  fh << var << " = np.array([\n";
80  for (std::ptrdiff_t r = 0; r < val.rows(); ++r)
81  {
82  if (r != 0)
83  {
84  fh << ",\n";
85  }
86  fh << " [";
87  for (std::ptrdiff_t c = 0; c < val.cols(); ++c)
88  {
89  if (c != 0)
90  {
91  fh << ", ";
92  }
93  fh << val(r, c);
94  }
95  fh << "]";
96  }
97  fh << "])\n";
98 }
99 
101 template <typename Derived>
102 void dump_py(std::ostream& fh, std::string const& var,
103  Eigen::ArrayBase<Derived> const& val)
104 {
105  dump_py(fh, var, val,
106  std::integral_constant<int, Derived::ColsAtCompileTime>{});
107 }
108 
110 template <typename Derived>
111 void dump_py(std::ostream& fh, std::string const& var,
112  Eigen::MatrixBase<Derived> const& val)
113 {
114  dump_py(fh, var, val.array());
115 }
116 
118 static const std::string msg_fatal =
119  "The local matrices M or K or the local vectors b assembled with the two "
120  "different Jacobian assemblers differ.";
121 
122 } // anonymous namespace
123 
124 namespace ProcessLib
125 {
127  LocalAssemblerInterface& local_assembler, double const t,
128  std::vector<double> const& local_x, std::vector<double> const& local_xdot,
129  const double dxdot_dx, const double dx_dx,
130  std::vector<double>& local_M_data, std::vector<double>& local_K_data,
131  std::vector<double>& local_b_data, std::vector<double>& local_Jac_data)
132 {
133  ++_counter;
134 
135  auto const num_dof = local_x.size();
136  auto to_mat = [num_dof](std::vector<double> const& data) {
137  if (data.empty())
138  {
139  return Eigen::Map<Eigen::Matrix<
140  double, Eigen::Dynamic, Eigen::Dynamic, Eigen::RowMajor> const>(
141  nullptr, 0, 0);
142  }
143 
144  return MathLib::toMatrix(data, num_dof, num_dof);
145  };
146 
147  // First assembly -- the one whose results will be added to the global
148  // equation system finally.
149  _asm1->assembleWithJacobian(local_assembler, t, local_x, local_xdot,
150  dxdot_dx, dx_dx, local_M_data, local_K_data,
151  local_b_data, local_Jac_data);
152 
153  auto const local_M1 = to_mat(local_M_data);
154  auto const local_K1 = to_mat(local_K_data);
155  auto const local_b1 = MathLib::toVector(local_b_data);
156 
157  std::vector<double> local_M_data2, local_K_data2, local_b_data2,
158  local_Jac_data2;
159 
160  // Second assembly -- used for checking only.
161  _asm2->assembleWithJacobian(local_assembler, t, local_x, local_xdot,
162  dxdot_dx, dx_dx, local_M_data2, local_K_data2,
163  local_b_data2, local_Jac_data2);
164 
165  auto const local_M2 = to_mat(local_M_data2);
166  auto const local_K2 = to_mat(local_K_data2);
167  auto const local_b2 = MathLib::toVector(local_b_data2);
168 
169  auto const local_Jac1 = MathLib::toMatrix(local_Jac_data, num_dof, num_dof);
170  auto const local_Jac2 =
171  MathLib::toMatrix(local_Jac_data2, num_dof, num_dof);
172 
173  auto const abs_diff = (local_Jac2 - local_Jac1).array().eval();
174  auto const rel_diff =
175  (abs_diff == 0.0)
176  .select(abs_diff,
177  2. * abs_diff /
178  (local_Jac1.cwiseAbs() + local_Jac2.cwiseAbs()).array())
179  .eval();
180 
181  auto const abs_diff_mask =
182  (abs_diff.abs() <= _abs_tol)
183  .select(decltype(abs_diff)::Zero(abs_diff.rows(), abs_diff.cols()),
184  decltype(abs_diff)::Ones(abs_diff.rows(), abs_diff.cols()))
185  .eval();
186  auto const rel_diff_mask =
187  (rel_diff.abs() <= _rel_tol)
188  .select(decltype(rel_diff)::Zero(rel_diff.rows(), rel_diff.cols()),
189  decltype(rel_diff)::Ones(rel_diff.rows(), rel_diff.cols()))
190  .eval();
191 
192  auto const abs_diff_OK = !abs_diff_mask.any();
193  auto const rel_diff_OK = !rel_diff_mask.any();
194 
195  std::ostringstream msg_tolerance;
196  bool tol_exceeded = true;
197  bool fatal_error = false;
198 
199  if (abs_diff_OK)
200  {
201  tol_exceeded = false;
202  }
203  else
204  {
205  msg_tolerance << "absolute tolerance of " << _abs_tol << " exceeded";
206  }
207 
208  if (rel_diff_OK)
209  {
210  tol_exceeded = false;
211  }
212  else
213  {
214  if (!msg_tolerance.str().empty())
215  {
216  msg_tolerance << " and ";
217  }
218 
219  msg_tolerance << "relative tolerance of " << _rel_tol << " exceeded";
220  }
221 
222  // basic consistency check if something went terribly wrong
223  auto check_equality = [&fatal_error](auto mat_or_vec1, auto mat_or_vec2) {
224  if (mat_or_vec1.size() == 0 || mat_or_vec2.size() == 0)
225  {
226  return;
227  }
228  if (mat_or_vec1.rows() != mat_or_vec2.rows() ||
229  mat_or_vec1.cols() != mat_or_vec2.cols())
230  {
231  fatal_error = true;
232  }
233  else if (((mat_or_vec1 - mat_or_vec2).array().cwiseAbs() >
234  std::numeric_limits<double>::epsilon())
235  .any())
236  {
237  fatal_error = true;
238  }
239  };
240 
241  check_equality(local_M1, local_M2);
242  check_equality(local_K1, local_K2);
243  check_equality(local_b1, local_b2);
244 
245  Eigen::VectorXd res1 = Eigen::VectorXd::Zero(num_dof);
246  if (local_M1.size() != 0)
247  {
248  res1.noalias() += local_M1 * MathLib::toVector(local_xdot);
249  }
250  if (local_K1.size() != 0)
251  {
252  res1.noalias() += local_K1 * MathLib::toVector(local_x);
253  }
254  if (local_b1.size() != 0)
255  {
256  res1.noalias() -= local_b1;
257  }
258 
259  Eigen::VectorXd res2 = Eigen::VectorXd::Zero(num_dof);
260  if (local_M2.size() != 0)
261  {
262  res2.noalias() += local_M2 * MathLib::toVector(local_xdot);
263  }
264  if (local_K2.size() != 0)
265  {
266  res2.noalias() += local_K2 * MathLib::toVector(local_x);
267  }
268  if (local_b2.size() != 0)
269  {
270  res2.noalias() -= local_b2;
271  }
272 
273  check_equality(res1, res2);
274 
275  if (tol_exceeded)
276  {
277  WARN("Compare Jacobians: %s", msg_tolerance.str().c_str());
278  }
279 
280  bool const output = tol_exceeded || fatal_error;
281 
282  if (output)
283  {
284  _log_file << "\n### counter: " << std::to_string(_counter)
285  << " (begin)\n";
286  }
287 
288  if (fatal_error)
289  {
290  _log_file << '\n'
291  << "#######################################################\n"
292  << "# FATAL ERROR: " << msg_fatal << '\n'
293  << "# You cannot expect any meaningful insights "
294  "from the Jacobian data printed below!\n"
295  << "# The reason for the mentioned differences "
296  "might be\n"
297  << "# (a) that the assembly routine has side "
298  "effects or\n"
299  << "# (b) that the assembly routines for M, K "
300  "and b themselves differ.\n"
301  << "#######################################################\n"
302  << '\n';
303  }
304 
305  if (tol_exceeded)
306  {
307  _log_file << "# " << msg_tolerance.str() << "\n\n";
308  }
309 
310  if (output)
311  {
312  dump_py(_log_file, "counter", _counter);
313  dump_py(_log_file, "num_dof", num_dof);
314  dump_py(_log_file, "abs_tol", _abs_tol);
315  dump_py(_log_file, "rel_tol", _rel_tol);
316 
317  _log_file << '\n';
318 
319  dump_py(_log_file, "local_x", local_x);
320  dump_py(_log_file, "local_x_dot", local_xdot);
321  dump_py(_log_file, "dxdot_dx", dxdot_dx);
322  dump_py(_log_file, "dx_dx", dx_dx);
323 
324  _log_file << '\n';
325 
326  dump_py(_log_file, "Jacobian_1", local_Jac1);
327  dump_py(_log_file, "Jacobian_2", local_Jac2);
328 
329  _log_file << '\n';
330 
331  _log_file << "# Jacobian_2 - Jacobian_1\n";
332  dump_py(_log_file, "abs_diff", abs_diff);
333  _log_file << "# Componentwise: 2 * abs_diff / (|Jacobian_1| + "
334  "|Jacobian_2|)\n";
335  dump_py(_log_file, "rel_diff", rel_diff);
336 
337  _log_file << '\n';
338 
339  _log_file << "# Masks: 0 ... tolerance met, 1 ... tolerance exceeded\n";
340  dump_py(_log_file, "abs_diff_mask", abs_diff_mask);
341  dump_py(_log_file, "rel_diff_mask", rel_diff_mask);
342 
343  _log_file << '\n';
344 
345  dump_py(_log_file, "M_1", local_M1);
346  dump_py(_log_file, "M_2", local_M2);
347  if (fatal_error && local_M1.size() == local_M2.size())
348  {
349  dump_py(_log_file, "delta_M", local_M2 - local_M1);
350  _log_file << '\n';
351  }
352 
353  dump_py(_log_file, "K_1", local_K1);
354  dump_py(_log_file, "K_2", local_K2);
355  if (fatal_error && local_K1.size() == local_K2.size())
356  {
357  dump_py(_log_file, "delta_K", local_K2 - local_K1);
358  _log_file << '\n';
359  }
360 
361  dump_py(_log_file, "b_1", local_b_data);
362  dump_py(_log_file, "b_2", local_b_data2);
363  if (fatal_error && local_b1.size() == local_b2.size())
364  {
365  dump_py(_log_file, "delta_b", local_b2 - local_b1);
366  _log_file << '\n';
367  }
368 
369  dump_py(_log_file, "res_1", res1);
370  dump_py(_log_file, "res_2", res2);
371  if (fatal_error)
372  {
373  dump_py(_log_file, "delta_res", res2 - res1);
374  }
375 
376  _log_file << '\n';
377 
378  _log_file << "### counter: " << std::to_string(_counter) << " (end)\n";
379  }
380 
381  if (fatal_error)
382  {
383  _log_file << std::flush;
384  OGS_FATAL("%s", msg_fatal.c_str());
385  }
386 
387  if (tol_exceeded && _fail_on_error)
388  {
389  _log_file << std::flush;
390  OGS_FATAL(
391  "OGS failed, because the two Jacobian implementations returned "
392  "different results.");
393  }
394 }
395 
396 std::unique_ptr<CompareJacobiansJacobianAssembler>
398 {
399  // TODO doc script corner case: Parameter could occur at different
400  // locations.
402  config.checkConfigParameter("type", "CompareJacobians");
403 
404  auto asm1 =
406  createJacobianAssembler(config.getConfigSubtree("jacobian_assembler"));
407 
408  auto asm2 = createJacobianAssembler(
410  config.getConfigSubtree("reference_jacobian_assembler"));
411 
413  auto const abs_tol = config.getConfigParameter<double>("abs_tol");
415  auto const rel_tol = config.getConfigParameter<double>("rel_tol");
416 
418  auto const fail_on_error = config.getConfigParameter<bool>("fail_on_error");
419 
421  auto const log_file = config.getConfigParameter<std::string>("log_file");
422 
423  return std::make_unique<CompareJacobiansJacobianAssembler>(
424  std::move(asm1), std::move(asm2), abs_tol, rel_tol, fail_on_error,
425  log_file);
426 }
427 
428 } // namespace ProcessLib
std::unique_ptr< AbstractJacobianAssembler > _asm1
void dump_py(std::ostream &fh, std::string const &var, Eigen::MatrixBase< Derived > const &val)
Dumps an Eigen matrix as a Python script snippet.
static const double r
std::unique_ptr< CompareJacobiansJacobianAssembler > createCompareJacobiansJacobianAssembler(BaseLib::ConfigTree const &config)
std::unique_ptr< AbstractJacobianAssembler > createJacobianAssembler(boost::optional< BaseLib::ConfigTree > const &config)
std::unique_ptr< AbstractJacobianAssembler > _asm2
Eigen::Map< const Vector > toVector(std::vector< double > const &data, Eigen::VectorXd::Index size)
Creates an Eigen mapped vector from the given data vector.
static const std::string msg_fatal
Will be printed if some consistency error is detected.
T getConfigParameter(std::string const &param) const
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
void checkConfigParameter(std::string const &param, T const &value) const
bool const _fail_on_error
Whether to abort if the tolerances are exceeded.
ConfigTree getConfigSubtree(std::string const &root) const
Definition: ConfigTree.cpp:150
void dump_py_vec(std::ostream &fh, std::string const &var, Vec const &val)
Dumps an arbitrary vector as a Python script snippet.
Eigen::Map< const Matrix > toMatrix(std::vector< double > const &data, Eigen::MatrixXd::Index rows, Eigen::MatrixXd::Index cols)
Definition: EigenMapTools.h:71
#define OGS_FATAL(fmt,...)
Definition: Error.h:63