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 31 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_x_prev, 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
 
std::unique_ptr< AbstractJacobianAssemblercopy () const override
 
- Public Member Functions inherited from ProcessLib::AbstractJacobianAssembler
virtual void assembleWithJacobianForStaggeredScheme (LocalAssemblerInterface &, double const, double const, Eigen::VectorXd const &, Eigen::VectorXd const &, 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.
 
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 34 of file CompareJacobiansJacobianAssembler.h.

41 : _asm1(std::move(asm1)),
42 _asm2(std::move(asm2)),
43 _abs_tol(abs_tol),
44 _rel_tol(rel_tol),
45 _fail_on_error(fail_on_error),
46 _log_file(log_file_path)
47 {
48 _log_file.precision(std::numeric_limits<double>::digits10);
49 _log_file << "#!/usr/bin/env python\n"
50 "import numpy as np\n"
51 "from numpy import nan\n"
52 << std::endl;
53 }
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_x_prev,
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 127 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_x_prev,
151 local_M_data, local_K_data, local_b_data,
152 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_x_prev,
165 local_M_data2, local_K_data2, local_b_data2,
166 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 auto const x = MathLib::toVector(local_x);
251 auto const x_dot = ((x - MathLib::toVector(local_x_prev)) / dt).eval();
252 if (local_M1.size() != 0)
253 {
254 res1.noalias() += local_M1 * x_dot;
255 }
256 if (local_K1.size() != 0)
257 {
258 res1.noalias() += local_K1 * x;
259 }
260 if (local_b1.size() != 0)
261 {
262 res1.noalias() -= local_b1;
263 }
264
265 Eigen::VectorXd res2 = Eigen::VectorXd::Zero(num_dof);
266 if (local_M2.size() != 0)
267 {
268 res2.noalias() += local_M2 * x_dot;
269 }
270 if (local_K2.size() != 0)
271 {
272 res2.noalias() += local_K2 * x;
273 }
274 if (local_b2.size() != 0)
275 {
276 res2.noalias() -= local_b2;
277 }
278
279 check_equality(res1, res2);
280
281 if (tol_exceeded)
282 {
283 WARN("Compare Jacobians: {:s}", msg_tolerance.str());
284 }
285
286 bool const output = tol_exceeded || fatal_error;
287
288 if (output)
289 {
290 _log_file << "\n### counter: " << std::to_string(_counter)
291 << " (begin)\n";
292 }
293
294 if (fatal_error)
295 {
296 _log_file << '\n'
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 "
302 "might be\n"
303 << "# (a) that the assembly routine has side "
304 "effects or\n"
305 << "# (b) that the assembly routines for M, K "
306 "and b themselves differ.\n"
307 << "#######################################################\n"
308 << '\n';
309 }
310
311 if (tol_exceeded)
312 {
313 _log_file << "# " << msg_tolerance.str() << "\n\n";
314 }
315
316 if (output)
317 {
318 dump_py(_log_file, "counter", _counter);
319 dump_py(_log_file, "num_dof", num_dof);
320 dump_py(_log_file, "abs_tol", _abs_tol);
321 dump_py(_log_file, "rel_tol", _rel_tol);
322
323 _log_file << '\n';
324
325 dump_py(_log_file, "local_x", local_x);
326 dump_py(_log_file, "local_x_prev", local_x_prev);
327 dump_py(_log_file, "dt", dt);
328
329 _log_file << '\n';
330
331 dump_py(_log_file, "Jacobian_1", local_Jac1);
332 dump_py(_log_file, "Jacobian_2", local_Jac2);
333
334 _log_file << '\n';
335
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| + "
339 "|Jacobian_2|)\n";
340 dump_py(_log_file, "rel_diff", rel_diff);
341
342 _log_file << '\n';
343
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);
347
348 _log_file << '\n';
349
350 dump_py(_log_file, "M_1", local_M1);
351 dump_py(_log_file, "M_2", local_M2);
352 if (fatal_error && local_M1.size() == local_M2.size())
353 {
354 dump_py(_log_file, "delta_M", local_M2 - local_M1);
355 _log_file << '\n';
356 }
357
358 dump_py(_log_file, "K_1", local_K1);
359 dump_py(_log_file, "K_2", local_K2);
360 if (fatal_error && local_K1.size() == local_K2.size())
361 {
362 dump_py(_log_file, "delta_K", local_K2 - local_K1);
363 _log_file << '\n';
364 }
365
366 dump_py(_log_file, "b_1", local_b_data);
367 dump_py(_log_file, "b_2", local_b_data2);
368 if (fatal_error && local_b1.size() == local_b2.size())
369 {
370 dump_py(_log_file, "delta_b", local_b2 - local_b1);
371 _log_file << '\n';
372 }
373
374 dump_py(_log_file, "res_1", res1);
375 dump_py(_log_file, "res_2", res2);
376 if (fatal_error)
377 {
378 dump_py(_log_file, "delta_res", res2 - res1);
379 }
380
381 _log_file << '\n';
382
383 _log_file << "### counter: " << std::to_string(_counter) << " (end)\n";
384 }
385
386 if (fatal_error)
387 {
388 _log_file << std::flush;
389 OGS_FATAL("{:s}", msg_fatal);
390 }
391
392 if (tol_exceeded && _fail_on_error)
393 {
394 _log_file << std::flush;
395 OGS_FATAL(
396 "OGS failed, because the two Jacobian implementations returned "
397 "different results.");
398 }
399}
#define OGS_FATAL(...)
Definition Error.h:26
void WARN(fmt::format_string< Args... > fmt, Args &&... args)
Definition Logging.h:40
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)
auto eval(Function &f, Tuples &... ts) -> typename detail::GetFunctionReturnType< decltype(&Function::eval)>::type
Definition Apply.h:274
void dump_py(std::ostream &fh, std::string const &var, double const val)
Dumps a double value 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, OGS_FATAL, MathLib::toMatrix(), MathLib::toVector(), and WARN().

◆ copy()

std::unique_ptr< AbstractJacobianAssembler > ProcessLib::CompareJacobiansJacobianAssembler::copy ( ) const
overridevirtual

Implements ProcessLib::AbstractJacobianAssembler.

Definition at line 402 of file CompareJacobiansJacobianAssembler.cpp.

403{
404 OGS_FATAL(
405 "CompareJacobiansJacobianAssembler should not be copied. This class "
406 "logs to a file, which would most certainly break after copying "
407 "(concurrent file access) with the current implementation.");
408}

References OGS_FATAL.

Member Data Documentation

◆ _abs_tol

double const ProcessLib::CompareJacobiansJacobianAssembler::_abs_tol
private

Definition at line 71 of file CompareJacobiansJacobianAssembler.h.

Referenced by assembleWithJacobian().

◆ _asm1

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

Definition at line 67 of file CompareJacobiansJacobianAssembler.h.

Referenced by assembleWithJacobian().

◆ _asm2

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

Definition at line 68 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 84 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 75 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 79 of file CompareJacobiansJacobianAssembler.h.

Referenced by CompareJacobiansJacobianAssembler(), and assembleWithJacobian().

◆ _rel_tol

double const ProcessLib::CompareJacobiansJacobianAssembler::_rel_tol
private

Definition at line 72 of file CompareJacobiansJacobianAssembler.h.

Referenced by assembleWithJacobian().


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