[PL] Added assembler for comparing Jacobians

ProcessLib/CompareJacobiansJacobianAssembler.cpp

+/**
+ * \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

ProcessLib/CompareJacobiansJacobianAssembler.h

+/**
+ * \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

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());
 }