21 void dump_py(std::ostream& fh, std::string
const& var,
double const val)
23 fh << var <<
" = " << val <<
'\n';
27 void dump_py(std::ostream& fh, std::string
const& var, std::size_t
const val)
29 fh << var <<
" = " << val <<
'\n';
33 template <
typename Vec>
34 void dump_py_vec(std::ostream& fh, std::string
const& var, Vec
const& val)
36 fh << var <<
" = np.array([";
37 for (decltype(val.size()) i = 0; i < val.size(); ++i)
58 void dump_py(std::ostream& fh, std::string
const& var,
59 std::vector<double>
const& val)
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> )
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> )
79 fh << var <<
" = np.array([\n";
80 for (std::ptrdiff_t
r = 0;
r < val.rows(); ++
r)
87 for (std::ptrdiff_t
c = 0;
c < val.cols(); ++
c)
101 template <
typename Derived>
102 void dump_py(std::ostream& fh, std::string
const& var,
103 Eigen::ArrayBase<Derived>
const& val)
106 std::integral_constant<int, Derived::ColsAtCompileTime>{});
110 template <
typename Derived>
111 void dump_py(std::ostream& fh, std::string
const& var,
112 Eigen::MatrixBase<Derived>
const& val)
119 "The local matrices M or K or the local vectors b assembled with the two "
120 "different Jacobian assemblers differ.";
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)
135 auto const num_dof = local_x.size();
136 auto to_mat = [num_dof](std::vector<double>
const& data)
140 return Eigen::Map<Eigen::Matrix<
141 double, Eigen::Dynamic, Eigen::Dynamic, Eigen::RowMajor>
const>(
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);
154 auto const local_M1 = to_mat(local_M_data);
155 auto const local_K1 = to_mat(local_K_data);
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;
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);
168 auto const local_M2 = to_mat(local_M_data2);
169 auto const local_K2 = to_mat(local_K_data2);
173 auto const local_Jac2 =
176 auto const abs_diff = (local_Jac2 - local_Jac1).array().eval();
177 auto const rel_diff =
181 (local_Jac1.cwiseAbs() + local_Jac2.cwiseAbs()).array())
184 auto const abs_diff_mask =
186 .select(decltype(abs_diff)::Zero(abs_diff.rows(), abs_diff.cols()),
187 decltype(abs_diff)::Ones(abs_diff.rows(), abs_diff.cols()))
189 auto const rel_diff_mask =
191 .select(decltype(rel_diff)::Zero(rel_diff.rows(), rel_diff.cols()),
192 decltype(rel_diff)::Ones(rel_diff.rows(), rel_diff.cols()))
195 auto const abs_diff_OK = !abs_diff_mask.any();
196 auto const rel_diff_OK = !rel_diff_mask.any();
198 std::ostringstream msg_tolerance;
199 bool tol_exceeded =
true;
200 bool fatal_error =
false;
204 tol_exceeded =
false;
208 msg_tolerance <<
"absolute tolerance of " <<
_abs_tol <<
" exceeded";
213 tol_exceeded =
false;
217 if (!msg_tolerance.str().empty())
219 msg_tolerance <<
" and ";
222 msg_tolerance <<
"relative tolerance of " <<
_rel_tol <<
" exceeded";
226 auto check_equality = [&fatal_error](
auto mat_or_vec1,
auto mat_or_vec2)
228 if (mat_or_vec1.size() == 0 || mat_or_vec2.size() == 0)
232 if (mat_or_vec1.rows() != mat_or_vec2.rows() ||
233 mat_or_vec1.cols() != mat_or_vec2.cols())
237 else if (((mat_or_vec1 - mat_or_vec2).array().cwiseAbs() >
238 std::numeric_limits<double>::epsilon())
245 check_equality(local_M1, local_M2);
246 check_equality(local_K1, local_K2);
247 check_equality(local_b1, local_b2);
249 Eigen::VectorXd res1 = Eigen::VectorXd::Zero(num_dof);
250 if (local_M1.size() != 0)
254 if (local_K1.size() != 0)
258 if (local_b1.size() != 0)
260 res1.noalias() -= local_b1;
263 Eigen::VectorXd res2 = Eigen::VectorXd::Zero(num_dof);
264 if (local_M2.size() != 0)
268 if (local_K2.size() != 0)
272 if (local_b2.size() != 0)
274 res2.noalias() -= local_b2;
277 check_equality(res1, res2);
281 WARN(
"Compare Jacobians: {:s}", msg_tolerance.str());
284 bool const output = tol_exceeded || fatal_error;
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 "
301 <<
"# (a) that the assembly routine has side "
303 <<
"# (b) that the assembly routines for M, K "
304 "and b themselves differ.\n"
305 <<
"#######################################################\n"
311 _log_file <<
"# " << msg_tolerance.str() <<
"\n\n";
335 _log_file <<
"# Jacobian_2 - Jacobian_1\n";
337 _log_file <<
"# Componentwise: 2 * abs_diff / (|Jacobian_1| + "
343 _log_file <<
"# Masks: 0 ... tolerance met, 1 ... tolerance exceeded\n";
351 if (fatal_error && local_M1.size() == local_M2.size())
359 if (fatal_error && local_K1.size() == local_K2.size())
367 if (fatal_error && local_b1.size() == local_b2.size())
395 "OGS failed, because the two Jacobian implementations returned "
396 "different results.");
400 std::unique_ptr<CompareJacobiansJacobianAssembler>
427 return std::make_unique<CompareJacobiansJacobianAssembler>(
428 std::move(asm1), std::move(asm2), abs_tol, rel_tol, fail_on_error,
void WARN(char const *fmt, Args const &... args)
void checkConfigParameter(std::string const ¶m, T const &value) const
T getConfigParameter(std::string const ¶m) const
ConfigTree getConfigSubtree(std::string const &root) const
bool const _fail_on_error
Whether to abort if the tolerances are exceeded.
std::unique_ptr< AbstractJacobianAssembler > _asm2
std::unique_ptr< AbstractJacobianAssembler > _asm1
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
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)
std::unique_ptr< CompareJacobiansJacobianAssembler > createCompareJacobiansJacobianAssembler(BaseLib::ConfigTree const &config)
std::unique_ptr< AbstractJacobianAssembler > createJacobianAssembler(std::optional< BaseLib::ConfigTree > const &config)
void dump_py_vec(std::ostream &fh, std::string const &var, Vec const &val)
Dumps an arbitrary vector as a Python script snippet.
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.