diff --git a/ProcessLib/CompareJacobiansJacobianAssembler.cpp b/ProcessLib/CompareJacobiansJacobianAssembler.cpp
new file mode 100644
index 0000000000000000000000000000000000000000..8c9972bd505039fdd042fb62482e3dfa6576ca3f
--- /dev/null
+++ b/ProcessLib/CompareJacobiansJacobianAssembler.cpp
@@ -0,0 +1,374 @@
+/**
+ * \copyright
+ * Copyright (c) 2012-2018, OpenGeoSys Community (http://www.opengeosys.org)
+ * Distributed under a Modified BSD License.
+ * See accompanying file LICENSE.txt or
+ * http://www.opengeosys.org/project/license
+ *
+ */
+
+#include "CompareJacobiansJacobianAssembler.h"
+
+#include <sstream>
+
+#include "MathLib/LinAlg/Eigen/EigenMapTools.h"
+
+#include "CreateJacobianAssembler.h"
+
+namespace
+{
+//! Dumps a \c double value as a Python script snippet.
+void dump_py(std::ofstream& fh, std::string const& var, double const val)
+{
+ fh << var << " = " << val << '\n';
+}
+
+//! Dumps a \c std::size_t value as a Python script snippet.
+void dump_py(std::ofstream& fh, std::string const& var, std::size_t const val)
+{
+ fh << var << " = " << val << '\n';
+}
+
+//! Dumps an Eigen vector as a Python script snippet.
+template <typename Vec>
+void dump_py_vec(std::ofstream& fh, std::string const& var, Vec const& val)
+{
+ fh << var << " = np.array([";
+ for (decltype(val.size()) i = 0; i < val.size(); ++i)
+ {
+ if (i != 0)
+ {
+ if (i % 8 == 0)
+ {
+ // Print at most eight entries on one line,
+ // indent with four spaces.
+ fh << ",\n ";
+ }
+ else
+ {
+ fh << ", ";
+ }
+ }
+ fh << val[i];
+ }
+ fh << "])\n";
+}
+
+//! Dumps a \c std::vector<double> as a Python script snippet.
+void dump_py(std::ofstream& fh, std::string const& var,
+ std::vector<double> const& val)
+{
+ dump_py_vec(fh, var, val);
+}
+
+template <typename Derived>
+void dump_py(std::ofstream& fh, std::string const& var,
+ Eigen::ArrayBase<Derived> const& val,
+ std::integral_constant<int, 1>)
+{
+ dump_py_vec(fh, var, val);
+}
+
+template <typename Derived, int ColsAtCompileTime>
+void dump_py(std::ofstream& fh, std::string const& var,
+ Eigen::ArrayBase<Derived> const& val,
+ std::integral_constant<int, ColsAtCompileTime>)
+{
+ fh << var << " = np.array([\n";
+ for (std::ptrdiff_t r = 0; r < val.rows(); ++r)
+ {
+ if (r != 0)
+ {
+ fh << ",\n";
+ }
+ fh << " [";
+ for (std::ptrdiff_t c = 0; c < val.cols(); ++c)
+ {
+ if (c != 0)
+ {
+ fh << ", ";
+ }
+ fh << val(r, c);
+ }
+ fh << "]";
+ }
+ fh << "])\n";
+}
+
+//! Dumps an Eigen array as a Python script snippet.
+template <typename Derived>
+void dump_py(std::ofstream& fh, std::string const& var,
+ Eigen::ArrayBase<Derived> const& val)
+{
+ dump_py(fh, var, val,
+ std::integral_constant<int, Derived::ColsAtCompileTime>{});
+}
+
+//! Dumps an Eigen matrix as a Python script snippet.
+template <typename Derived>
+void dump_py(std::ofstream& fh, std::string const& var,
+ Eigen::MatrixBase<Derived> const& val)
+{
+ dump_py(fh, var, val.array());
+}
+
+//! Will be printed if some consistency error is detected.
+static const std::string msg_fatal =
+ "The local matrices M or K or the local vectors b assembled with the two "
+ "different Jacobian assemblers differ.";
+
+} // anonymous namespace
+
+namespace ProcessLib
+{
+void CompareJacobiansJacobianAssembler::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)
+{
+ ++_counter;
+
+ auto const num_dof = local_x.size();
+ auto to_mat = [num_dof](std::vector<double> const& data) {
+ if (data.empty())
+ {
+ return Eigen::Map<Eigen::Matrix<
+ double, Eigen::Dynamic, Eigen::Dynamic, Eigen::RowMajor> const>(
+ nullptr, 0, 0);
+ }
+
+ return MathLib::toMatrix(data, num_dof, num_dof);
+ };
+
+ // First assembly -- the one whose results will be added to the global
+ // equation system finally.
+ _asm1->assembleWithJacobian(local_assembler, t, local_x, local_xdot,
+ dxdot_dx, dx_dx, local_M_data, local_K_data,
+ local_b_data, local_Jac_data);
+
+ auto const local_M1 = to_mat(local_M_data);
+ auto const local_K1 = to_mat(local_K_data);
+ auto const local_b1 = MathLib::toVector(local_b_data);
+
+ std::vector<double> local_M_data2, local_K_data2, local_b_data2,
+ local_Jac_data2;
+
+ // Second assembly -- used for checking only.
+ _asm2->assembleWithJacobian(local_assembler, t, local_x, local_xdot,
+ dxdot_dx, dx_dx, local_M_data2, local_K_data2,
+ local_b_data2, local_Jac_data2);
+
+ auto const local_M2 = to_mat(local_M_data2);
+ auto const local_K2 = to_mat(local_K_data2);
+ auto const local_b2 = MathLib::toVector(local_b_data2);
+
+ auto const local_Jac1 = MathLib::toMatrix(local_Jac_data, num_dof, num_dof);
+ auto const local_Jac2 =
+ MathLib::toMatrix(local_Jac_data2, num_dof, num_dof);
+
+ auto const abs_diff = (local_Jac2 - local_Jac1).array().eval();
+ auto const rel_diff =
+ (abs_diff == 0.0)
+ .select(abs_diff,
+ 2. * abs_diff /
+ (local_Jac1.cwiseAbs() + local_Jac2.cwiseAbs()).array())
+ .eval();
+
+ // the !(... <= ...) construct handles NaN correctly.
+ auto const abs_diff_mask = (!(abs_diff.cwiseAbs() <= _abs_tol)).eval();
+ auto const rel_diff_mask = (!(rel_diff.cwiseAbs() <= _rel_tol)).eval();
+
+ auto const abs_diff_OK = !abs_diff_mask.any();
+ auto const rel_diff_OK = !rel_diff_mask.any();
+
+ std::ostringstream msg_tolerance;
+ bool tol_exceeded = true;
+ bool fatal_error = false;
+
+ if (abs_diff_OK)
+ {
+ tol_exceeded = false;
+ }
+ else
+ {
+ msg_tolerance << "absolute tolerance of " << _abs_tol << " exceeded";
+ }
+
+ if (rel_diff_OK)
+ {
+ tol_exceeded = false;
+ }
+ else
+ {
+ if (!msg_tolerance.str().empty())
+ {
+ msg_tolerance << " and ";
+ }
+
+ msg_tolerance << "relative tolerance of " << _rel_tol << " exceeded";
+ }
+
+ // basic consistency check if something went terribly wrong
+ auto check_equality = [&fatal_error](auto mat_or_vec1, auto mat_or_vec2) {
+ if (mat_or_vec1.rows() != mat_or_vec2.rows() ||
+ mat_or_vec1.cols() != mat_or_vec2.cols())
+ {
+ fatal_error = true;
+ }
+ else if (((mat_or_vec1 - mat_or_vec2).array().cwiseAbs() >
+ std::numeric_limits<double>::epsilon())
+ .any())
+ {
+ fatal_error = true;
+ }
+ };
+ check_equality(local_M1, local_M2);
+ check_equality(local_K1, local_K2);
+ check_equality(local_b1, local_b2);
+
+ if (tol_exceeded)
+ {
+ WARN("Compare Jacobians: %s", msg_tolerance.str().c_str());
+ }
+
+ bool const output = tol_exceeded || fatal_error;
+
+ if (output)
+ {
+ _log_file << "\n### counter: " << std::to_string(_counter)
+ << " (begin)\n";
+ }
+
+ if (fatal_error)
+ {
+ _log_file << '\n'
+ << "#######################################################\n"
+ << "# FATAL ERROR: " << msg_fatal << '\n'
+ << "# You cannot expect any meaningful insights "
+ "from the Jacobian data printed below!\n"
+ << "# The reason for the mentioned differences "
+ "might be\n"
+ << "# (a) that the assembly routine has side "
+ "effects or\n"
+ << "# (b) that the assembly routines for M, K "
+ "and b themselves differ.\n"
+ << "#######################################################\n"
+ << '\n';
+ }
+
+ if (tol_exceeded)
+ {
+ _log_file << "# " << msg_tolerance.str() << "\n\n";
+ }
+
+ if (output)
+ {
+ dump_py(_log_file, "counter", _counter);
+ dump_py(_log_file, "num_dof", num_dof);
+ dump_py(_log_file, "abs_tol", _abs_tol);
+ dump_py(_log_file, "rel_tol", _rel_tol);
+
+ _log_file << '\n';
+
+ dump_py(_log_file, "local_x", local_x);
+ dump_py(_log_file, "local_x_dot", local_xdot);
+ dump_py(_log_file, "dxdot_dx", dxdot_dx);
+ dump_py(_log_file, "dx_dx", dx_dx);
+
+ _log_file << '\n';
+
+ dump_py(_log_file, "Jacobian_1", local_Jac1);
+ dump_py(_log_file, "Jacobian_2", local_Jac2);
+
+ _log_file << '\n';
+
+ _log_file << "# Jacobian_2 - Jacobian_1\n";
+ dump_py(_log_file, "abs_diff", abs_diff);
+ _log_file << "# componentwise: 2 * abs_diff / (|Jacobian_1| + "
+ "|Jacobian_2|)\n";
+ dump_py(_log_file, "rel_diff", rel_diff);
+
+ _log_file << '\n';
+
+ _log_file << "# Masks: 0 ... tolerance met, 1 ... tolerance exceeded\n";
+ dump_py(_log_file, "abs_diff_mask", abs_diff_mask);
+ dump_py(_log_file, "rel_diff_mask", rel_diff_mask);
+
+ _log_file << '\n';
+
+ dump_py(_log_file, "M_1", local_M1);
+ dump_py(_log_file, "M_2", local_M2);
+ if (fatal_error)
+ {
+ dump_py(_log_file, "delta_M", local_M2 - local_M1);
+ _log_file << '\n';
+ }
+
+ dump_py(_log_file, "K_1", local_K1);
+ dump_py(_log_file, "K_2", local_K2);
+ if (fatal_error)
+ {
+ dump_py(_log_file, "delta_K", local_K2 - local_K1);
+ _log_file << '\n';
+ }
+
+ dump_py(_log_file, "b_1", local_b_data);
+ dump_py(_log_file, "b_2", local_b_data2);
+ if (fatal_error)
+ {
+ dump_py(_log_file, "delta_b", local_b2 - local_b1);
+ }
+
+ _log_file << '\n';
+
+ _log_file << "### counter: " << std::to_string(_counter) << " (end)\n";
+ }
+
+ if (fatal_error)
+ {
+ OGS_FATAL("%s", msg_fatal.c_str());
+ }
+
+ if (tol_exceeded && _fail_on_error)
+ {
+ OGS_FATAL(
+ "OGS failed, because the two Jacobian implementations returned "
+ "different results.");
+ }
+}
+
+std::unique_ptr<CompareJacobiansJacobianAssembler>
+createCompareJacobiansJacobianAssembler(BaseLib::ConfigTree const& config)
+{
+ // TODO doc script corner case: Parameter could occur at different
+ // locations.
+ //! \ogs_file_param{prj__processes__process__jacobian_assembler__type}
+ config.checkConfigParameter("type", "CompareJacobians");
+
+ //! \ogs_file_param{prj__processes__process__jacobian_assembler__CompareJacobians__jacobian_assembler}
+ auto asm1 =
+ createJacobianAssembler(config.getConfigSubtree("jacobian_assembler"));
+
+ //! \ogs_file_param{prj__processes__process__jacobian_assembler__CompareJacobians__reference_jacobian_assembler}
+ auto asm2 = createJacobianAssembler(
+ config.getConfigSubtree("reference_jacobian_assembler"));
+
+ //! \ogs_file_param{prj__processes__process__jacobian_assembler__CompareJacobians__abs_tol}
+ auto const abs_tol = config.getConfigParameter<double>("abs_tol");
+ //! \ogs_file_param{prj__processes__process__jacobian_assembler__CompareJacobians__rel_tol}
+ auto const rel_tol = config.getConfigParameter<double>("rel_tol");
+
+ //! \ogs_file_param{prj__processes__process__jacobian_assembler__CompareJacobians__fail_on_error}
+ auto const fail_on_error = config.getConfigParameter<bool>("fail_on_error");
+
+ //! \ogs_file_param{prj__processes__process__jacobian_assembler__CompareJacobians__log_file}
+ auto const log_file = config.getConfigParameter<std::string>("log_file");
+
+ return std::make_unique<CompareJacobiansJacobianAssembler>(
+ std::move(asm1), std::move(asm2), abs_tol, rel_tol, fail_on_error,
+ log_file);
+}
+
+} // namespace ProcessLib
diff --git a/ProcessLib/CompareJacobiansJacobianAssembler.h b/ProcessLib/CompareJacobiansJacobianAssembler.h
new file mode 100644
index 0000000000000000000000000000000000000000..71e752afa881f9a53a8593b4811c4f8c2b3084fa
--- /dev/null
+++ b/ProcessLib/CompareJacobiansJacobianAssembler.h
@@ -0,0 +1,87 @@
+/**
+ * \copyright
+ * Copyright (c) 2012-2018, OpenGeoSys Community (http://www.opengeosys.org)
+ * Distributed under a Modified BSD License.
+ * See accompanying file LICENSE.txt or
+ * http://www.opengeosys.org/project/license
+ *
+ */
+
+#pragma once
+
+#include <fstream>
+#include <limits>
+#include <memory>
+#include "AbstractJacobianAssembler.h"
+
+namespace BaseLib
+{
+class ConfigTree;
+} // namespace BaseLib
+
+namespace ProcessLib
+{
+//! 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.
+class CompareJacobiansJacobianAssembler final : public AbstractJacobianAssembler
+{
+public:
+ 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)
+ : _asm1(std::move(asm1)),
+ _asm2(std::move(asm2)),
+ _abs_tol(abs_tol),
+ _rel_tol(rel_tol),
+ _fail_on_error(fail_on_error),
+ _log_file(log_file_path)
+ {
+ _log_file.precision(std::numeric_limits<double>::digits10);
+ _log_file << "#!/usr/bin/env python\n"
+ "import numpy as np\n"
+ "from numpy import nan\n"
+ << std::endl;
+ }
+
+ 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;
+
+private:
+ std::unique_ptr<AbstractJacobianAssembler> _asm1;
+ std::unique_ptr<AbstractJacobianAssembler> _asm2;
+
+ // TODO change to matrix
+ double const _abs_tol;
+ double const _rel_tol;
+
+ //! Whether to abort if the tolerances are exceeded.
+ bool const _fail_on_error;
+
+ //! Path where a Python script will be placed, which contains information
+ //! about exceeded tolerances and assembled local matrices.
+ std::ofstream _log_file;
+
+ //! Counter used for identifying blocks in the \c _log_file. Is be
+ //! incremented upon each call of the assembly routine, i.e., for each
+ //! element in each iteration etc.
+ std::size_t _counter = 0;
+};
+
+std::unique_ptr<CompareJacobiansJacobianAssembler>
+createCompareJacobiansJacobianAssembler(BaseLib::ConfigTree const& config);
+
+} // namespace ProcessLib
diff --git a/ProcessLib/CreateJacobianAssembler.cpp b/ProcessLib/CreateJacobianAssembler.cpp
index d677daed7d0f9cea3e236f08633bc3cecd910686..39f80d66313126ff8c95291034986b2f43c77e0b 100644
--- a/ProcessLib/CreateJacobianAssembler.cpp
+++ b/ProcessLib/CreateJacobianAssembler.cpp
@@ -12,6 +12,7 @@
#include "AnalyticalJacobianAssembler.h"
#include "CentralDifferencesJacobianAssembler.h"
+#include "CompareJacobiansJacobianAssembler.h"
namespace ProcessLib
{
@@ -32,6 +33,10 @@ std::unique_ptr<AbstractJacobianAssembler> createJacobianAssembler(
{
return createCentralDifferencesJacobianAssembler(*config);
}
+ if (type == "CompareJacobians")
+ {
+ return createCompareJacobiansJacobianAssembler(*config);
+ }
OGS_FATAL("Unknown Jacobian assembler type: `%s'.", type.c_str());
}