129 std::vector<double>
const& local_x, std::vector<double>
const& local_x_prev,
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_x_prev,
151 local_M_data, local_K_data, local_b_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_x_prev,
165 local_M_data2, local_K_data2, local_b_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);
252 if (local_M1.size() != 0)
254 res1.noalias() += local_M1 * x_dot;
256 if (local_K1.size() != 0)
258 res1.noalias() += local_K1 * x;
260 if (local_b1.size() != 0)
262 res1.noalias() -= local_b1;
265 Eigen::VectorXd res2 = Eigen::VectorXd::Zero(num_dof);
266 if (local_M2.size() != 0)
268 res2.noalias() += local_M2 * x_dot;
270 if (local_K2.size() != 0)
272 res2.noalias() += local_K2 * x;
274 if (local_b2.size() != 0)
276 res2.noalias() -= local_b2;
279 check_equality(res1, res2);
283 WARN(
"Compare Jacobians: {:s}", msg_tolerance.str());
286 bool const output = tol_exceeded || fatal_error;
297 <<
"#######################################################\n"
298 <<
"# FATAL ERROR: " << msg_fatal <<
'\n'
299 <<
"# You cannot expect any meaningful insights "
300 "from the Jacobian data printed below!\n"
301 <<
"# The reason for the mentioned differences "
303 <<
"# (a) that the assembly routine has side "
305 <<
"# (b) that the assembly routines for M, K "
306 "and b themselves differ.\n"
307 <<
"#######################################################\n"
313 _log_file <<
"# " << msg_tolerance.str() <<
"\n\n";
326 dump_py(
_log_file,
"local_x_prev", local_x_prev);
331 dump_py(
_log_file,
"Jacobian_1", local_Jac1);
332 dump_py(
_log_file,
"Jacobian_2", local_Jac2);
336 _log_file <<
"# Jacobian_2 - Jacobian_1\n";
337 dump_py(
_log_file,
"abs_diff", abs_diff);
338 _log_file <<
"# Componentwise: 2 * abs_diff / (|Jacobian_1| + "
340 dump_py(
_log_file,
"rel_diff", rel_diff);
344 _log_file <<
"# Masks: 0 ... tolerance met, 1 ... tolerance exceeded\n";
345 dump_py(
_log_file,
"abs_diff_mask", abs_diff_mask);
346 dump_py(
_log_file,
"rel_diff_mask", rel_diff_mask);
352 if (fatal_error && local_M1.size() == local_M2.size())
354 dump_py(
_log_file,
"delta_M", local_M2 - local_M1);
360 if (fatal_error && local_K1.size() == local_K2.size())
362 dump_py(
_log_file,
"delta_K", local_K2 - local_K1);
367 dump_py(
_log_file,
"b_2", local_b_data2);
368 if (fatal_error && local_b1.size() == local_b2.size())
370 dump_py(
_log_file,
"delta_b", local_b2 - local_b1);
378 dump_py(
_log_file,
"delta_res", res2 - res1);
396 "OGS failed, because the two Jacobian implementations returned "
397 "different results.");