ProcessLib::CompareJacobiansJacobianAssembler Class Reference

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:

Collaboration diagram for ProcessLib::CompareJacobiansJacobianAssembler:

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< AbstractJacobianAssembler >	copy () 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.

        : _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;
    }

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.

{
    ++_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, dt, local_x, local_x_prev,
                                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;
    std::vector<double> local_K_data2;
    std::vector<double> local_b_data2;
    std::vector<double> local_Jac_data2;
 
    // Second assembly -- used for checking only.
    _asm2->assembleWithJacobian(local_assembler, t, dt, local_x, local_x_prev,
                                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();
 
    auto const abs_diff_mask =
        (abs_diff.abs() <= _abs_tol)
            .select(decltype(abs_diff)::Zero(abs_diff.rows(), abs_diff.cols()),
                    decltype(abs_diff)::Ones(abs_diff.rows(), abs_diff.cols()))
            .eval();
    auto const rel_diff_mask =
        (rel_diff.abs() <= _rel_tol)
            .select(decltype(rel_diff)::Zero(rel_diff.rows(), rel_diff.cols()),
                    decltype(rel_diff)::Ones(rel_diff.rows(), rel_diff.cols()))
            .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.size() == 0 || mat_or_vec2.size() == 0)
        {
            return;
        }
        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);
 
    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_M1.size() != 0)
    {
        res1.noalias() += local_M1 * x_dot;
    }
    if (local_K1.size() != 0)
    {
        res1.noalias() += local_K1 * x;
    }
    if (local_b1.size() != 0)
    {
        res1.noalias() -= local_b1;
    }
 
    Eigen::VectorXd res2 = Eigen::VectorXd::Zero(num_dof);
    if (local_M2.size() != 0)
    {
        res2.noalias() += local_M2 * x_dot;
    }
    if (local_K2.size() != 0)
    {
        res2.noalias() += local_K2 * x;
    }
    if (local_b2.size() != 0)
    {
        res2.noalias() -= local_b2;
    }
 
    check_equality(res1, res2);
 
    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)
                  << " (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_prev", local_x_prev);
        dump_py(_log_file, "dt", dt);
 
        _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 && local_M1.size() == local_M2.size())
        {
            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 && local_K1.size() == local_K2.size())
        {
            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 && 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.");
    }
}

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.

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

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:

• ProcessLib/CompareJacobiansJacobianAssembler.h
• ProcessLib/CompareJacobiansJacobianAssembler.cpp