OGS 6.2.1-97-g73d1aeda3
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> /*unused*/)
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> /*unused*/)
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 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;
158  std::vector<double> local_K_data2;
159  std::vector<double> local_b_data2;
160  std::vector<double> local_Jac_data2;
161 
162  // Second assembly -- used for checking only.
163  _asm2->assembleWithJacobian(local_assembler, t, local_x, local_xdot,
164  dxdot_dx, dx_dx, local_M_data2, local_K_data2,
165  local_b_data2, local_Jac_data2);
166 
167  auto const local_M2 = to_mat(local_M_data2);
168  auto const local_K2 = to_mat(local_K_data2);
169  auto const local_b2 = MathLib::toVector(local_b_data2);
170 
171  auto const local_Jac1 = MathLib::toMatrix(local_Jac_data, num_dof, num_dof);
172  auto const local_Jac2 =
173  MathLib::toMatrix(local_Jac_data2, num_dof, num_dof);
174 
175  auto const abs_diff = (local_Jac2 - local_Jac1).array().eval();
176  auto const rel_diff =
177  (abs_diff == 0.0)
178  .select(abs_diff,
179  2. * abs_diff /
180  (local_Jac1.cwiseAbs() + local_Jac2.cwiseAbs()).array())
181  .eval();
182 
183  auto const abs_diff_mask =
184  (abs_diff.abs() <= _abs_tol)
185  .select(decltype(abs_diff)::Zero(abs_diff.rows(), abs_diff.cols()),
186  decltype(abs_diff)::Ones(abs_diff.rows(), abs_diff.cols()))
187  .eval();
188  auto const rel_diff_mask =
189  (rel_diff.abs() <= _rel_tol)
190  .select(decltype(rel_diff)::Zero(rel_diff.rows(), rel_diff.cols()),
191  decltype(rel_diff)::Ones(rel_diff.rows(), rel_diff.cols()))
192  .eval();
193 
194  auto const abs_diff_OK = !abs_diff_mask.any();
195  auto const rel_diff_OK = !rel_diff_mask.any();
196 
197  std::ostringstream msg_tolerance;
198  bool tol_exceeded = true;
199  bool fatal_error = false;
200 
201  if (abs_diff_OK)
202  {
203  tol_exceeded = false;
204  }
205  else
206  {
207  msg_tolerance << "absolute tolerance of " << _abs_tol << " exceeded";
208  }
209 
210  if (rel_diff_OK)
211  {
212  tol_exceeded = false;
213  }
214  else
215  {
216  if (!msg_tolerance.str().empty())
217  {
218  msg_tolerance << " and ";
219  }
220 
221  msg_tolerance << "relative tolerance of " << _rel_tol << " exceeded";
222  }
223 
224  // basic consistency check if something went terribly wrong
225  auto check_equality = [&fatal_error](auto mat_or_vec1, auto mat_or_vec2) {
226  if (mat_or_vec1.size() == 0 || mat_or_vec2.size() == 0)
227  {
228  return;
229  }
230  if (mat_or_vec1.rows() != mat_or_vec2.rows() ||
231  mat_or_vec1.cols() != mat_or_vec2.cols())
232  {
233  fatal_error = true;
234  }
235  else if (((mat_or_vec1 - mat_or_vec2).array().cwiseAbs() >
236  std::numeric_limits<double>::epsilon())
237  .any())
238  {
239  fatal_error = true;
240  }
241  };
242 
243  check_equality(local_M1, local_M2);
244  check_equality(local_K1, local_K2);
245  check_equality(local_b1, local_b2);
246 
247  Eigen::VectorXd res1 = Eigen::VectorXd::Zero(num_dof);
248  if (local_M1.size() != 0)
249  {
250  res1.noalias() += local_M1 * MathLib::toVector(local_xdot);
251  }
252  if (local_K1.size() != 0)
253  {
254  res1.noalias() += local_K1 * MathLib::toVector(local_x);
255  }
256  if (local_b1.size() != 0)
257  {
258  res1.noalias() -= local_b1;
259  }
260 
261  Eigen::VectorXd res2 = Eigen::VectorXd::Zero(num_dof);
262  if (local_M2.size() != 0)
263  {
264  res2.noalias() += local_M2 * MathLib::toVector(local_xdot);
265  }
266  if (local_K2.size() != 0)
267  {
268  res2.noalias() += local_K2 * MathLib::toVector(local_x);
269  }
270  if (local_b2.size() != 0)
271  {
272  res2.noalias() -= local_b2;
273  }
274 
275  check_equality(res1, res2);
276 
277  if (tol_exceeded)
278  {
279  WARN("Compare Jacobians: %s", msg_tolerance.str().c_str());
280  }
281 
282  bool const output = tol_exceeded || fatal_error;
283 
284  if (output)
285  {
286  _log_file << "\n### counter: " << std::to_string(_counter)
287  << " (begin)\n";
288  }
289 
290  if (fatal_error)
291  {
292  _log_file << '\n'
293  << "#######################################################\n"
294  << "# FATAL ERROR: " << msg_fatal << '\n'
295  << "# You cannot expect any meaningful insights "
296  "from the Jacobian data printed below!\n"
297  << "# The reason for the mentioned differences "
298  "might be\n"
299  << "# (a) that the assembly routine has side "
300  "effects or\n"
301  << "# (b) that the assembly routines for M, K "
302  "and b themselves differ.\n"
303  << "#######################################################\n"
304  << '\n';
305  }
306 
307  if (tol_exceeded)
308  {
309  _log_file << "# " << msg_tolerance.str() << "\n\n";
310  }
311 
312  if (output)
313  {
314  dump_py(_log_file, "counter", _counter);
315  dump_py(_log_file, "num_dof", num_dof);
316  dump_py(_log_file, "abs_tol", _abs_tol);
317  dump_py(_log_file, "rel_tol", _rel_tol);
318 
319  _log_file << '\n';
320 
321  dump_py(_log_file, "local_x", local_x);
322  dump_py(_log_file, "local_x_dot", local_xdot);
323  dump_py(_log_file, "dxdot_dx", dxdot_dx);
324  dump_py(_log_file, "dx_dx", dx_dx);
325 
326  _log_file << '\n';
327 
328  dump_py(_log_file, "Jacobian_1", local_Jac1);
329  dump_py(_log_file, "Jacobian_2", local_Jac2);
330 
331  _log_file << '\n';
332 
333  _log_file << "# Jacobian_2 - Jacobian_1\n";
334  dump_py(_log_file, "abs_diff", abs_diff);
335  _log_file << "# Componentwise: 2 * abs_diff / (|Jacobian_1| + "
336  "|Jacobian_2|)\n";
337  dump_py(_log_file, "rel_diff", rel_diff);
338 
339  _log_file << '\n';
340 
341  _log_file << "# Masks: 0 ... tolerance met, 1 ... tolerance exceeded\n";
342  dump_py(_log_file, "abs_diff_mask", abs_diff_mask);
343  dump_py(_log_file, "rel_diff_mask", rel_diff_mask);
344 
345  _log_file << '\n';
346 
347  dump_py(_log_file, "M_1", local_M1);
348  dump_py(_log_file, "M_2", local_M2);
349  if (fatal_error && local_M1.size() == local_M2.size())
350  {
351  dump_py(_log_file, "delta_M", local_M2 - local_M1);
352  _log_file << '\n';
353  }
354 
355  dump_py(_log_file, "K_1", local_K1);
356  dump_py(_log_file, "K_2", local_K2);
357  if (fatal_error && local_K1.size() == local_K2.size())
358  {
359  dump_py(_log_file, "delta_K", local_K2 - local_K1);
360  _log_file << '\n';
361  }
362 
363  dump_py(_log_file, "b_1", local_b_data);
364  dump_py(_log_file, "b_2", local_b_data2);
365  if (fatal_error && local_b1.size() == local_b2.size())
366  {
367  dump_py(_log_file, "delta_b", local_b2 - local_b1);
368  _log_file << '\n';
369  }
370 
371  dump_py(_log_file, "res_1", res1);
372  dump_py(_log_file, "res_2", res2);
373  if (fatal_error)
374  {
375  dump_py(_log_file, "delta_res", res2 - res1);
376  }
377 
378  _log_file << '\n';
379 
380  _log_file << "### counter: " << std::to_string(_counter) << " (end)\n";
381  }
382 
383  if (fatal_error)
384  {
385  _log_file << std::flush;
386  OGS_FATAL("%s", msg_fatal.c_str());
387  }
388 
389  if (tol_exceeded && _fail_on_error)
390  {
391  _log_file << std::flush;
392  OGS_FATAL(
393  "OGS failed, because the two Jacobian implementations returned "
394  "different results.");
395  }
396 }
397 
398 std::unique_ptr<CompareJacobiansJacobianAssembler>
400 {
401  // TODO doc script corner case: Parameter could occur at different
402  // locations.
404  config.checkConfigParameter("type", "CompareJacobians");
405 
406  auto asm1 =
408  createJacobianAssembler(config.getConfigSubtree("jacobian_assembler"));
409 
410  auto asm2 = createJacobianAssembler(
412  config.getConfigSubtree("reference_jacobian_assembler"));
413 
415  auto const abs_tol = config.getConfigParameter<double>("abs_tol");
417  auto const rel_tol = config.getConfigParameter<double>("rel_tol");
418 
420  auto const fail_on_error = config.getConfigParameter<bool>("fail_on_error");
421 
423  auto const log_file = config.getConfigParameter<std::string>("log_file");
424 
425  return std::make_unique<CompareJacobiansJacobianAssembler>(
426  std::move(asm1), std::move(asm2), abs_tol, rel_tol, fail_on_error,
427  log_file);
428 }
429 
430 } // 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.
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