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 24 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_b_data, std::vector< double > &local_Jac_data) override
std::unique_ptr< AbstractJacobianAssembler >	copy () const override
Public Member Functions inherited from ProcessLib::AbstractJacobianAssembler
	AbstractJacobianAssembler (std::vector< double > const &&absolute_epsilons)
	AbstractJacobianAssembler ()
virtual void	assembleWithJacobianForStaggeredScheme (LocalAssemblerInterface &, double const, double const, Eigen::VectorXd const &, Eigen::VectorXd const &, int const, std::vector< double > &, std::vector< double > &)
virtual	~AbstractJacobianAssembler ()=default
void	checkPerturbationSize (int const max_non_deformation_dofs_per_node) const
void	setNonDeformationComponentIDs (std::vector< int > const &non_deformation_component_ids)
void	setNonDeformationComponentIDsNoSizeCheck (std::vector< int > const &non_deformation_component_ids)
auto	getVariableComponentEpsilonsView () const
auto	isPerturbationEnabled () const

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

Additional Inherited Members
Protected Attributes inherited from ProcessLib::AbstractJacobianAssembler
std::vector< int >	non_deformation_component_ids_
std::vector< double > const	absolute_epsilons_

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 27 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>::max_digits10);
        _log_file << "#!/usr/bin/env python\n"
                     "import numpy as np\n"
                     "from numpy import nan\n"
                  << std::endl;
    }

References _abs_tol, _asm1, _asm2, _fail_on_error, _log_file, and _rel_tol.

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_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 120 of file CompareJacobiansJacobianAssembler.cpp.

{
    ++_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(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(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 = (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_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_b1.size() != 0)
    {
        res1.noalias() -= local_b1;
    }
 
    Eigen::VectorXd res2 = Eigen::VectorXd::Zero(num_dof);
    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 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, "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 341 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 62 of file CompareJacobiansJacobianAssembler.h.

Referenced by CompareJacobiansJacobianAssembler(), and assembleWithJacobian().

_asm1

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

private

Definition at line 58 of file CompareJacobiansJacobianAssembler.h.

Referenced by CompareJacobiansJacobianAssembler(), and assembleWithJacobian().

_asm2

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

private

Definition at line 59 of file CompareJacobiansJacobianAssembler.h.

Referenced by CompareJacobiansJacobianAssembler(), and 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 75 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 66 of file CompareJacobiansJacobianAssembler.h.

Referenced by CompareJacobiansJacobianAssembler(), and 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 70 of file CompareJacobiansJacobianAssembler.h.

Referenced by CompareJacobiansJacobianAssembler(), and assembleWithJacobian().

_rel_tol

double const ProcessLib::CompareJacobiansJacobianAssembler::_rel_tol

private

Definition at line 63 of file CompareJacobiansJacobianAssembler.h.

Referenced by CompareJacobiansJacobianAssembler(), and assembleWithJacobian().

---

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

• CompareJacobiansJacobianAssembler.h
• CompareJacobiansJacobianAssembler.cpp