OGS
ProcessLib::CompareJacobiansJacobianAssembler Class Referencefinal

Detailed Description

Assembles the Jacobian matrix using two different Jacobian assemblers and compares the assembled local Jacobian matrices.

If the provided tolerances are exceeded, debugging information is logged in the form of a Python script.

Definition at line 30 of file CompareJacobiansJacobianAssembler.h.

#include <CompareJacobiansJacobianAssembler.h>

Inheritance diagram for ProcessLib::CompareJacobiansJacobianAssembler:
[legend]
Collaboration diagram for ProcessLib::CompareJacobiansJacobianAssembler:
[legend]

Public Member Functions

 CompareJacobiansJacobianAssembler (std::unique_ptr< AbstractJacobianAssembler > &&asm1, std::unique_ptr< AbstractJacobianAssembler > &&asm2, double abs_tol, double rel_tol, bool fail_on_error, std::string const &log_file_path)
 
void assembleWithJacobian (LocalAssemblerInterface &local_assembler, double const t, double const dt, 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
 
- Public Member Functions inherited from ProcessLib::AbstractJacobianAssembler
virtual void assembleWithJacobianForStaggeredScheme (LocalAssemblerInterface &, double const, double const, Eigen::VectorXd const &, Eigen::VectorXd const &, const double, const double, int const, std::vector< double > &, std::vector< double > &, std::vector< double > &, std::vector< double > &)
 
virtual ~AbstractJacobianAssembler ()=default
 

Private Attributes

std::unique_ptr< AbstractJacobianAssembler_asm1
 
std::unique_ptr< AbstractJacobianAssembler_asm2
 
double const _abs_tol
 
double const _rel_tol
 
bool const _fail_on_error
 Whether to abort if the tolerances are exceeded. More...
 
std::ofstream _log_file
 
std::size_t _counter = 0
 

Constructor & Destructor Documentation

◆ CompareJacobiansJacobianAssembler()

ProcessLib::CompareJacobiansJacobianAssembler::CompareJacobiansJacobianAssembler ( std::unique_ptr< AbstractJacobianAssembler > &&  asm1,
std::unique_ptr< AbstractJacobianAssembler > &&  asm2,
double  abs_tol,
double  rel_tol,
bool  fail_on_error,
std::string const &  log_file_path 
)
inline

Definition at line 33 of file CompareJacobiansJacobianAssembler.h.

40  : _asm1(std::move(asm1)),
41  _asm2(std::move(asm2)),
42  _abs_tol(abs_tol),
43  _rel_tol(rel_tol),
44  _fail_on_error(fail_on_error),
45  _log_file(log_file_path)
46  {
47  _log_file.precision(std::numeric_limits<double>::digits10);
48  _log_file << "#!/usr/bin/env python\n"
49  "import numpy as np\n"
50  "from numpy import nan\n"
51  << std::endl;
52  }
bool const _fail_on_error
Whether to abort if the tolerances are exceeded.
std::unique_ptr< AbstractJacobianAssembler > _asm2
std::unique_ptr< AbstractJacobianAssembler > _asm1

References _log_file.

Member Function Documentation

◆ assembleWithJacobian()

void ProcessLib::CompareJacobiansJacobianAssembler::assembleWithJacobian ( LocalAssemblerInterface local_assembler,
double const  t,
double const  dt,
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 
)
overridevirtual

Assembles the Jacobian, the matrices \(M\) and \(K\), and the vector \(b\).

Implements ProcessLib::AbstractJacobianAssembler.

Definition at line 126 of file CompareJacobiansJacobianAssembler.cpp.

132 {
133  ++_counter;
134 
135  auto const num_dof = local_x.size();
136  auto to_mat = [num_dof](std::vector<double> const& data)
137  {
138  if (data.empty())
139  {
140  return Eigen::Map<Eigen::Matrix<
141  double, Eigen::Dynamic, Eigen::Dynamic, Eigen::RowMajor> const>(
142  nullptr, 0, 0);
143  }
144 
145  return MathLib::toMatrix(data, num_dof, num_dof);
146  };
147 
148  // First assembly -- the one whose results will be added to the global
149  // equation system finally.
150  _asm1->assembleWithJacobian(local_assembler, t, dt, local_x, local_xdot,
151  dxdot_dx, dx_dx, local_M_data, local_K_data,
152  local_b_data, local_Jac_data);
153 
154  auto const local_M1 = to_mat(local_M_data);
155  auto const local_K1 = to_mat(local_K_data);
156  auto const local_b1 = MathLib::toVector(local_b_data);
157 
158  std::vector<double> local_M_data2;
159  std::vector<double> local_K_data2;
160  std::vector<double> local_b_data2;
161  std::vector<double> local_Jac_data2;
162 
163  // Second assembly -- used for checking only.
164  _asm2->assembleWithJacobian(local_assembler, t, dt, local_x, local_xdot,
165  dxdot_dx, dx_dx, local_M_data2, local_K_data2,
166  local_b_data2, local_Jac_data2);
167 
168  auto const local_M2 = to_mat(local_M_data2);
169  auto const local_K2 = to_mat(local_K_data2);
170  auto const local_b2 = MathLib::toVector(local_b_data2);
171 
172  auto const local_Jac1 = MathLib::toMatrix(local_Jac_data, num_dof, num_dof);
173  auto const local_Jac2 =
174  MathLib::toMatrix(local_Jac_data2, num_dof, num_dof);
175 
176  auto const abs_diff = (local_Jac2 - local_Jac1).array().eval();
177  auto const rel_diff =
178  (abs_diff == 0.0)
179  .select(abs_diff,
180  2. * abs_diff /
181  (local_Jac1.cwiseAbs() + local_Jac2.cwiseAbs()).array())
182  .eval();
183 
184  auto const abs_diff_mask =
185  (abs_diff.abs() <= _abs_tol)
186  .select(decltype(abs_diff)::Zero(abs_diff.rows(), abs_diff.cols()),
187  decltype(abs_diff)::Ones(abs_diff.rows(), abs_diff.cols()))
188  .eval();
189  auto const rel_diff_mask =
190  (rel_diff.abs() <= _rel_tol)
191  .select(decltype(rel_diff)::Zero(rel_diff.rows(), rel_diff.cols()),
192  decltype(rel_diff)::Ones(rel_diff.rows(), rel_diff.cols()))
193  .eval();
194 
195  auto const abs_diff_OK = !abs_diff_mask.any();
196  auto const rel_diff_OK = !rel_diff_mask.any();
197 
198  std::ostringstream msg_tolerance;
199  bool tol_exceeded = true;
200  bool fatal_error = false;
201 
202  if (abs_diff_OK)
203  {
204  tol_exceeded = false;
205  }
206  else
207  {
208  msg_tolerance << "absolute tolerance of " << _abs_tol << " exceeded";
209  }
210 
211  if (rel_diff_OK)
212  {
213  tol_exceeded = false;
214  }
215  else
216  {
217  if (!msg_tolerance.str().empty())
218  {
219  msg_tolerance << " and ";
220  }
221 
222  msg_tolerance << "relative tolerance of " << _rel_tol << " exceeded";
223  }
224 
225  // basic consistency check if something went terribly wrong
226  auto check_equality = [&fatal_error](auto mat_or_vec1, auto mat_or_vec2)
227  {
228  if (mat_or_vec1.size() == 0 || mat_or_vec2.size() == 0)
229  {
230  return;
231  }
232  if (mat_or_vec1.rows() != mat_or_vec2.rows() ||
233  mat_or_vec1.cols() != mat_or_vec2.cols())
234  {
235  fatal_error = true;
236  }
237  else if (((mat_or_vec1 - mat_or_vec2).array().cwiseAbs() >
238  std::numeric_limits<double>::epsilon())
239  .any())
240  {
241  fatal_error = true;
242  }
243  };
244 
245  check_equality(local_M1, local_M2);
246  check_equality(local_K1, local_K2);
247  check_equality(local_b1, local_b2);
248 
249  Eigen::VectorXd res1 = Eigen::VectorXd::Zero(num_dof);
250  if (local_M1.size() != 0)
251  {
252  res1.noalias() += local_M1 * MathLib::toVector(local_xdot);
253  }
254  if (local_K1.size() != 0)
255  {
256  res1.noalias() += local_K1 * MathLib::toVector(local_x);
257  }
258  if (local_b1.size() != 0)
259  {
260  res1.noalias() -= local_b1;
261  }
262 
263  Eigen::VectorXd res2 = Eigen::VectorXd::Zero(num_dof);
264  if (local_M2.size() != 0)
265  {
266  res2.noalias() += local_M2 * MathLib::toVector(local_xdot);
267  }
268  if (local_K2.size() != 0)
269  {
270  res2.noalias() += local_K2 * MathLib::toVector(local_x);
271  }
272  if (local_b2.size() != 0)
273  {
274  res2.noalias() -= local_b2;
275  }
276 
277  check_equality(res1, res2);
278 
279  if (tol_exceeded)
280  {
281  WARN("Compare Jacobians: {:s}", msg_tolerance.str());
282  }
283 
284  bool const output = tol_exceeded || fatal_error;
285 
286  if (output)
287  {
288  _log_file << "\n### counter: " << std::to_string(_counter)
289  << " (begin)\n";
290  }
291 
292  if (fatal_error)
293  {
294  _log_file << '\n'
295  << "#######################################################\n"
296  << "# FATAL ERROR: " << msg_fatal << '\n'
297  << "# You cannot expect any meaningful insights "
298  "from the Jacobian data printed below!\n"
299  << "# The reason for the mentioned differences "
300  "might be\n"
301  << "# (a) that the assembly routine has side "
302  "effects or\n"
303  << "# (b) that the assembly routines for M, K "
304  "and b themselves differ.\n"
305  << "#######################################################\n"
306  << '\n';
307  }
308 
309  if (tol_exceeded)
310  {
311  _log_file << "# " << msg_tolerance.str() << "\n\n";
312  }
313 
314  if (output)
315  {
316  dump_py(_log_file, "counter", _counter);
317  dump_py(_log_file, "num_dof", num_dof);
318  dump_py(_log_file, "abs_tol", _abs_tol);
319  dump_py(_log_file, "rel_tol", _rel_tol);
320 
321  _log_file << '\n';
322 
323  dump_py(_log_file, "local_x", local_x);
324  dump_py(_log_file, "local_x_dot", local_xdot);
325  dump_py(_log_file, "dxdot_dx", dxdot_dx);
326  dump_py(_log_file, "dx_dx", dx_dx);
327 
328  _log_file << '\n';
329 
330  dump_py(_log_file, "Jacobian_1", local_Jac1);
331  dump_py(_log_file, "Jacobian_2", local_Jac2);
332 
333  _log_file << '\n';
334 
335  _log_file << "# Jacobian_2 - Jacobian_1\n";
336  dump_py(_log_file, "abs_diff", abs_diff);
337  _log_file << "# Componentwise: 2 * abs_diff / (|Jacobian_1| + "
338  "|Jacobian_2|)\n";
339  dump_py(_log_file, "rel_diff", rel_diff);
340 
341  _log_file << '\n';
342 
343  _log_file << "# Masks: 0 ... tolerance met, 1 ... tolerance exceeded\n";
344  dump_py(_log_file, "abs_diff_mask", abs_diff_mask);
345  dump_py(_log_file, "rel_diff_mask", rel_diff_mask);
346 
347  _log_file << '\n';
348 
349  dump_py(_log_file, "M_1", local_M1);
350  dump_py(_log_file, "M_2", local_M2);
351  if (fatal_error && local_M1.size() == local_M2.size())
352  {
353  dump_py(_log_file, "delta_M", local_M2 - local_M1);
354  _log_file << '\n';
355  }
356 
357  dump_py(_log_file, "K_1", local_K1);
358  dump_py(_log_file, "K_2", local_K2);
359  if (fatal_error && local_K1.size() == local_K2.size())
360  {
361  dump_py(_log_file, "delta_K", local_K2 - local_K1);
362  _log_file << '\n';
363  }
364 
365  dump_py(_log_file, "b_1", local_b_data);
366  dump_py(_log_file, "b_2", local_b_data2);
367  if (fatal_error && local_b1.size() == local_b2.size())
368  {
369  dump_py(_log_file, "delta_b", local_b2 - local_b1);
370  _log_file << '\n';
371  }
372 
373  dump_py(_log_file, "res_1", res1);
374  dump_py(_log_file, "res_2", res2);
375  if (fatal_error)
376  {
377  dump_py(_log_file, "delta_res", res2 - res1);
378  }
379 
380  _log_file << '\n';
381 
382  _log_file << "### counter: " << std::to_string(_counter) << " (end)\n";
383  }
384 
385  if (fatal_error)
386  {
387  _log_file << std::flush;
388  OGS_FATAL("{:s}", msg_fatal);
389  }
390 
391  if (tol_exceeded && _fail_on_error)
392  {
393  _log_file << std::flush;
394  OGS_FATAL(
395  "OGS failed, because the two Jacobian implementations returned "
396  "different results.");
397  }
398 }
#define OGS_FATAL(...)
Definition: Error.h:26
void WARN(char const *fmt, Args const &... args)
Definition: Logging.h:37
Eigen::Map< const Vector > toVector(std::vector< double > const &data, Eigen::VectorXd::Index size)
Creates an Eigen mapped vector from the given data vector.
Eigen::Map< const Matrix > toMatrix(std::vector< double > const &data, Eigen::MatrixXd::Index rows, Eigen::MatrixXd::Index cols)
Definition: EigenMapTools.h:72
void dump_py(std::ostream &fh, std::string const &var, Eigen::MatrixBase< Derived > const &val)
Dumps an Eigen matrix as a Python script snippet.
const std::string msg_fatal
Will be printed if some consistency error is detected.

References _abs_tol, _asm1, _asm2, _counter, _fail_on_error, _log_file, _rel_tol, anonymous_namespace{CompareJacobiansJacobianAssembler.cpp}::dump_py(), anonymous_namespace{CompareJacobiansJacobianAssembler.cpp}::msg_fatal, OGS_FATAL, MathLib::toMatrix(), MathLib::toVector(), and WARN().

Member Data Documentation

◆ _abs_tol

double const ProcessLib::CompareJacobiansJacobianAssembler::_abs_tol
private

Definition at line 69 of file CompareJacobiansJacobianAssembler.h.

Referenced by assembleWithJacobian().

◆ _asm1

std::unique_ptr<AbstractJacobianAssembler> ProcessLib::CompareJacobiansJacobianAssembler::_asm1
private

Definition at line 65 of file CompareJacobiansJacobianAssembler.h.

Referenced by assembleWithJacobian().

◆ _asm2

std::unique_ptr<AbstractJacobianAssembler> ProcessLib::CompareJacobiansJacobianAssembler::_asm2
private

Definition at line 66 of file CompareJacobiansJacobianAssembler.h.

Referenced by assembleWithJacobian().

◆ _counter

std::size_t ProcessLib::CompareJacobiansJacobianAssembler::_counter = 0
private

Counter used for identifying blocks in the _log_file. It is incremented upon each call of the assembly routine, i.e., for each element in each iteration etc.

Definition at line 82 of file CompareJacobiansJacobianAssembler.h.

Referenced by assembleWithJacobian().

◆ _fail_on_error

bool const ProcessLib::CompareJacobiansJacobianAssembler::_fail_on_error
private

Whether to abort if the tolerances are exceeded.

Definition at line 73 of file CompareJacobiansJacobianAssembler.h.

Referenced by assembleWithJacobian().

◆ _log_file

std::ofstream ProcessLib::CompareJacobiansJacobianAssembler::_log_file
private

Path where a Python script will be placed, which contains information about exceeded tolerances and assembled local matrices.

Definition at line 77 of file CompareJacobiansJacobianAssembler.h.

Referenced by CompareJacobiansJacobianAssembler(), and assembleWithJacobian().

◆ _rel_tol

double const ProcessLib::CompareJacobiansJacobianAssembler::_rel_tol
private

Definition at line 70 of file CompareJacobiansJacobianAssembler.h.

Referenced by assembleWithJacobian().


The documentation for this class was generated from the following files: