CompareJacobiansJacobianAssembler.cpp

Go to the documentation of this file.

// SPDX-FileCopyrightText: Copyright (c) OpenGeoSys Community (opengeosys.org)
// SPDX-License-Identifier: BSD-3-Clause
 
#include "CompareJacobiansJacobianAssembler.h"
 
#include <fstream>
#include <limits>
#include <sstream>
 
#ifdef _OPENMP
#include <omp.h>
#endif
 
#include "BaseLib/ConfigTree.h"
#include "CreateJacobianAssembler.h"
#include "MathLib/LinAlg/Eigen/EigenMapTools.h"
 
namespace
{
template <typename T>
    requires std::integral<T> || std::floating_point<T>
void dump_py(std::ostream& fh, std::string const& var, T const val)
{
    fh << var << " = " << val << '\n';
}

template <typename Vec>
void dump_py_vec(std::ostream& 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";
}

void dump_py(std::ostream& fh, std::string const& var,
             std::vector<double> const& val)
{
    dump_py_vec(fh, var, val);
}

template <typename Derived>
void dump_py(std::ostream& fh, std::string const& var,
             Eigen::ArrayBase<Derived> const& val,
             std::integral_constant<int, 1> /*unused*/)
{
    dump_py_vec(fh, var, val);
}

template <typename Derived, int ColsAtCompileTime>
void dump_py(std::ostream& fh, std::string const& var,
             Eigen::ArrayBase<Derived> const& val,
             std::integral_constant<int, ColsAtCompileTime> /*unused*/)
{
    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";
}

template <typename Derived>
void dump_py(std::ostream& fh, std::string const& var,
             Eigen::ArrayBase<Derived> const& val)
{
    dump_py(fh, var, val,
            std::integral_constant<int, Derived::ColsAtCompileTime>{});
}

template <typename Derived>
void dump_py(std::ostream& fh, std::string const& var,
             Eigen::MatrixBase<Derived> const& val)
{
    dump_py(fh, var, val.array());
}

const std::string msg_fatal =
    "The local matrices M or K or the local vectors b assembled with the two "
    "different Jacobian assemblers differ.";

auto signedAbsDiffRelDiff(auto const& mat_or_vec1, auto const& mat_or_vec2)
{
    auto abs_diff = (mat_or_vec2 - mat_or_vec1).array().eval();
    auto const rel_diff =
        (abs_diff == 0.0)
            .select(
                abs_diff,
                2. * abs_diff /
                    (mat_or_vec1.cwiseAbs() + mat_or_vec2.cwiseAbs()).array())
            .eval();
    return std::pair{std::move(abs_diff), std::move(rel_diff)};
}
 
template <typename MatOrVec>
auto oneIfAboveThresholdElseZero(MatOrVec&& mat_or_vec, double const threshold)
{
    return (mat_or_vec <= threshold)
        .select(MatOrVec::Zero(mat_or_vec.rows(), mat_or_vec.cols()),
                MatOrVec::Ones(mat_or_vec.rows(), mat_or_vec.cols()))
        .eval();
}
 
// basic consistency check if something went terribly wrong
bool isSimilar(auto const& mat_or_vec1, auto const& mat_or_vec2,
               double const abs_tol, double const rel_tol)
{
    if (mat_or_vec1.size() == 0 || mat_or_vec2.size() == 0)
    {
        // either size is 0, ignore
        return true;
    }
    if (mat_or_vec1.rows() != mat_or_vec2.rows() ||
        mat_or_vec1.cols() != mat_or_vec2.cols())
    {
        return false;
    }
 
    auto const [abs_diff, rel_diff] =
        signedAbsDiffRelDiff(mat_or_vec1, mat_or_vec2);
    auto const abs_tol_exceeded = abs_diff.abs() > abs_tol;
    auto const rel_tol_exceeded = rel_diff.abs() > rel_tol;
 
    // similar if for no entry abs and rel tols are exceeded at the same time
    return !(abs_tol_exceeded && rel_tol_exceeded).any();
}
 
}  // anonymous namespace
 
namespace ProcessLib
{
namespace detail
{
struct CompareJacobiansJacobianAssemblerImpl
{
    friend class ProcessLib::CompareJacobiansJacobianAssembler;
 
    CompareJacobiansJacobianAssemblerImpl(
        std::unique_ptr<AbstractJacobianAssembler>&& asm1,
        std::unique_ptr<AbstractJacobianAssembler>&& asm2,
        double abs_tol_Jac,
        double rel_tol_Jac,
        double abs_tol_res,
        double rel_tol_res,
        bool fail_on_error,
        std::string const& log_file_path)
        : asm1_{std::move(asm1)},
          asm2_{std::move(asm2)},
          abs_tol_Jac_{abs_tol_Jac},
          rel_tol_Jac_{rel_tol_Jac},
          abs_tol_res_{abs_tol_res},
          rel_tol_res_{rel_tol_res},
          fail_on_error_{fail_on_error},
          log_file_{log_file_path}
    {
        log_file_.precision(std::numeric_limits<double>::max_digits10);
        log_file_ << "#!/usr/bin/env python\n"
                     "import numpy as np\n"
                     "from numpy import nan\n"
                  << std::endl;
    }
 
    void assembleWithJacobian(std::size_t const mesh_item_id,
                              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_b_data,
                              std::vector<double>& local_Jac_data);
 
private:
    std::unique_ptr<AbstractJacobianAssembler> asm1_;
    std::unique_ptr<AbstractJacobianAssembler> asm2_;
 
    double const abs_tol_Jac_;
    double const rel_tol_Jac_;
    double const abs_tol_res_;
    double const rel_tol_res_;

    bool const fail_on_error_;

    std::ofstream log_file_;

    std::ptrdiff_t counter_ = -1;
 
    unsigned iter_ = 0;
};
 
void CompareJacobiansJacobianAssemblerImpl::assembleWithJacobian(
    std::size_t const mesh_item_id, 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_b_data,
    std::vector<double>& local_Jac_data)
{
    ++counter_;
 
    auto const num_dof = local_x.size();
 
    // First assembly -- the one whose results will be added to the global
    // equation system finally.
    asm1_->assembleWithJacobian(mesh_item_id, local_assembler, t, dt, local_x,
                                local_x_prev, local_b_data, local_Jac_data);
 
    auto const local_b1 = MathLib::toVector(local_b_data);
 
    std::vector<double> local_b_data2;
    std::vector<double> local_Jac_data2;
 
    // Second assembly -- used for checking only.
    asm2_->assembleWithJacobian(mesh_item_id, local_assembler, t, dt, local_x,
                                local_x_prev, local_b_data2, local_Jac_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, rel_diff] =
        signedAbsDiffRelDiff(local_Jac1, local_Jac2);
 
    auto const abs_diff_mask =
        oneIfAboveThresholdElseZero(abs_diff.abs(), abs_tol_Jac_);
    auto const rel_diff_mask =
        oneIfAboveThresholdElseZero(rel_diff.abs(), rel_tol_Jac_);
 
    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_Jac_
                      << " exceeded";
    }
 
    if (rel_diff_OK)
    {
        tol_exceeded = false;
    }
    else
    {
        if (!msg_tolerance.str().empty())
        {
            msg_tolerance << " and ";
        }
 
        msg_tolerance << "relative tolerance of " << rel_tol_Jac_
                      << " exceeded";
    }
 
    fatal_error |= !isSimilar(local_b1, local_b2, abs_tol_res_, rel_tol_res_);
 
    Eigen::VectorXd res1 = Eigen::VectorXd::Zero(num_dof);
    auto const x = MathLib::toVector(local_x);
    auto const x_dot = ((x - MathLib::toVector(local_x_prev)) / dt).eval();
    if (local_b1.size() != 0)
    {
        res1.noalias() -= local_b1;
    }
 
    Eigen::VectorXd res2 = Eigen::VectorXd::Zero(num_dof);
    if (local_b2.size() != 0)
    {
        res2.noalias() -= local_b2;
    }
 
    fatal_error |= !isSimilar(res1, res2, abs_tol_res_, rel_tol_res_);
 
    if (tol_exceeded)
    {
        WARN("Compare Jacobians: {:s}", msg_tolerance.str());
    }
 
    bool const output = tol_exceeded || fatal_error;
 
    if (output)
    {
        log_file_ << "\n### counter: " << std::to_string(counter_)
                  << ", t: " << t << ", element_id: " << mesh_item_id
                  << " (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 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_, "nonlinear_iteration", iter_);
        dump_py(log_file_, "t", t);
        dump_py(log_file_, "dt", dt);
        dump_py(log_file_, "element_id", mesh_item_id);
 
        log_file_ << '\n';
 
        dump_py(log_file_, "num_dof", num_dof);
        dump_py(log_file_, "abs_tol", abs_tol_Jac_);
        dump_py(log_file_, "rel_tol", rel_tol_Jac_);
 
        log_file_ << '\n';
 
        dump_py(log_file_, "local_x", local_x);
        dump_py(log_file_, "local_x_prev", local_x_prev);
 
        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_ << "# max(|abs_diff|) = " << abs_diff.abs().maxCoeff()
                  << '\n';
        log_file_ << "# Componentwise: 2 * abs_diff / (|Jacobian_1| + "
                     "|Jacobian_2|)\n";
        dump_py(log_file_, "rel_diff", rel_diff);
        log_file_ << "# max(|rel_diff|) = " << rel_diff.abs().maxCoeff()
                  << '\n';
 
        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_, "b_1", local_b_data);
        dump_py(log_file_, "b_2", local_b_data2);
        if (fatal_error && local_b1.size() == local_b2.size())
        {
            dump_py(log_file_, "delta_b", local_b2 - local_b1);
            log_file_ << '\n';
        }
 
        dump_py(log_file_, "res_1", res1);
        dump_py(log_file_, "res_2", res2);
        if (fatal_error)
        {
            dump_py(log_file_, "delta_res", res2 - res1);
        }
 
        log_file_ << '\n';
 
        log_file_ << "### counter: " << std::to_string(counter_) << " (end)\n";
    }
 
    if (fatal_error)
    {
        log_file_ << std::flush;
        OGS_FATAL("{:s}", msg_fatal);
    }
 
    if (tol_exceeded && fail_on_error_)
    {
        log_file_ << std::flush;
        OGS_FATAL(
            "OGS failed, because the two Jacobian implementations returned "
            "different results.");
    }
}
}  // namespace detail
 
CompareJacobiansJacobianAssembler::CompareJacobiansJacobianAssembler(
    std::unique_ptr<AbstractJacobianAssembler>&& asm1,
    std::unique_ptr<AbstractJacobianAssembler>&& asm2, double abs_tol_Jac,
    double rel_tol_Jac, double abs_tol_res, double rel_tol_res,
    bool fail_on_error, std::string const& log_file_path)
    : impl_{std::make_shared<detail::CompareJacobiansJacobianAssemblerImpl>(
          std::move(asm1), std::move(asm2), abs_tol_Jac, rel_tol_Jac,
          abs_tol_res, rel_tol_res, fail_on_error, log_file_path)}
{
}
 
CompareJacobiansJacobianAssembler::CompareJacobiansJacobianAssembler(
    std::shared_ptr<detail::CompareJacobiansJacobianAssemblerImpl> impl,
    CompareJacobiansJacobianAssembler::Key /*key*/)
    : impl_{std::move(impl)}
{
}
 
void CompareJacobiansJacobianAssembler::assembleWithJacobian(
    std::size_t const mesh_item_id, 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_b_data,
    std::vector<double>& local_Jac_data)
{
    impl_->assembleWithJacobian(mesh_item_id, local_assembler, t, dt, local_x,
                                local_x_prev, local_b_data, local_Jac_data);
}
 
std::unique_ptr<AbstractJacobianAssembler>
CompareJacobiansJacobianAssembler::copy() const
{
#ifdef _OPENMP
    if (omp_get_thread_num() != 0)
    {
        OGS_FATAL(
            "CompareJacobiansJacobianAssembler cannot be used concurrently. "
            "Please restrict yourself to one assembly thread "
            "(OGS_ASM_THREADS=1).");
    }
#endif
 
    return std::make_unique<CompareJacobiansJacobianAssembler>(impl_, Key{});
}
 
void CompareJacobiansJacobianAssembler::checkPerturbationSize(
    int const max_non_deformation_dofs_per_node) const
{
    impl_->asm1_->checkPerturbationSize(max_non_deformation_dofs_per_node);
    impl_->asm2_->checkPerturbationSize(max_non_deformation_dofs_per_node);
}
 
void CompareJacobiansJacobianAssembler::setNonDeformationComponentIDs(
    std::vector<int> const& non_deformation_component_ids)
{
    impl_->asm1_->setNonDeformationComponentIDs(non_deformation_component_ids);
    impl_->asm2_->setNonDeformationComponentIDs(non_deformation_component_ids);
}
 
void CompareJacobiansJacobianAssembler::
    setNonDeformationComponentIDsNoSizeCheck(
        std::vector<int> const& non_deformation_component_ids)
{
    impl_->asm1_->setNonDeformationComponentIDsNoSizeCheck(
        non_deformation_component_ids);
    impl_->asm2_->setNonDeformationComponentIDsNoSizeCheck(
        non_deformation_component_ids);
}
 
bool CompareJacobiansJacobianAssembler::needsPicardAssembly() const
{
    return impl_->asm1_->needsPicardAssembly() ||
           impl_->asm2_->needsPicardAssembly();
}
 
void CompareJacobiansJacobianAssembler::preIteration(const unsigned iter)
{
    impl_->iter_ = iter;
}
 
std::unique_ptr<CompareJacobiansJacobianAssembler>
createCompareJacobiansJacobianAssembler(BaseLib::ConfigTree const& config)
{
    // TODO doc script corner case: Parameter could occur at different
    // locations.
    config.checkConfigParameter("type", "CompareJacobians");
 
    auto asm1 =
        createJacobianAssembler(config.getConfigSubtree("jacobian_assembler"));
 
    auto asm2 = createJacobianAssembler(
        config.getConfigSubtree("reference_jacobian_assembler"));

    auto const abs_tol = config.getConfigParameter<double>("abs_tol");
    auto const rel_tol = config.getConfigParameter<double>("rel_tol");

    auto const abs_tol_res = config.getConfigParameter<double>("abs_tol_res");
    auto const rel_tol_res = config.getConfigParameter<double>("rel_tol_res");

    auto const fail_on_error = config.getConfigParameter<bool>("fail_on_error");

    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, abs_tol_res,
        rel_tol_res, fail_on_error, log_file);
}
}  // namespace ProcessLib