Skip to content

GitLab

Commit 58de6df4 authored by GILLES Sebastien's avatar GILLES Sebastien
Browse files

#1335 Remove entirely the macro MOREFEM_CHECK_NAN_AND_INF; the few tests... 

#1335 Remove entirely the macro MOREFEM_CHECK_NAN_AND_INF; the few tests checked were actually too constraining for some models (hyperelastic models may lead to infinite in some Newton steps, in which case the Newton diverges ad another attempt might be made).

I put in Utilities/ThirdParty new functions to check whether a floating point, a Seldon vector or a Seldon matrix includes some inf or nan; it is not currently in use but a unit test has been added to ensure their correct behaviour.
parent 24d881bd
Hide whitespace changes
Inline Side-by-side
Showing with 181 additions and 157 deletions
CMakeLists.txt
View file @ 58de6df4
......@@ -20,11 +20,6 @@ if(MOREFEM_EXTENDED_TIME_KEEP)
add_definitions(-DMOREFEM_EXTENDED_TIME_KEEP)
endif()
if(MOREFEM_CHECK_NAN_AND_INF)
message("ADDING FLAG MOREFEM_CHECK_NAN_AND_INF")
add_definitions(-DMOREFEM_CHECK_NAN_AND_INF)
endif()
set_property(GLOBAL PROPERTY USE_FOLDERS ON)
set(CMAKE_POSITION_INDEPENDENT_CODE ON)
......
MoReFEM.xcodeproj/project.pbxproj
View file @ 58de6df4
......@@ -4792,6 +4792,9 @@
BEE734D01F75432D001D01A9 /* demo_input_data.lua */ = {isa = PBXFileReference; fileEncoding = 4; lastKnownFileType = text; path = demo_input_data.lua; sourceTree = "<group>"; };
BEE734D11F75432D001D01A9 /* InputParameterList.hpp */ = {isa = PBXFileReference; fileEncoding = 4; lastKnownFileType = sourcecode.cpp.h; path = InputParameterList.hpp; sourceTree = "<group>"; };
BEE734D51F7550F7001D01A9 /* test.cpp */ = {isa = PBXFileReference; fileEncoding = 4; lastKnownFileType = sourcecode.cpp.cpp; path = test.cpp; sourceTree = "<group>"; };
BEE79C6A2136F43000388FED /* CMakeLists.txt */ = {isa = PBXFileReference; lastKnownFileType = text; path = CMakeLists.txt; sourceTree = "<group>"; };
BEE79C6B2136F43000388FED /* test.cpp */ = {isa = PBXFileReference; lastKnownFileType = sourcecode.cpp.cpp; path = test.cpp; sourceTree = "<group>"; };
BEE79C6C2136F51E00388FED /* CMakeLists.txt */ = {isa = PBXFileReference; fileEncoding = 4; lastKnownFileType = text; path = CMakeLists.txt; sourceTree = "<group>"; };
BEE8A2821C90488700CD25F0 /* Viscosity.cpp */ = {isa = PBXFileReference; fileEncoding = 4; lastKnownFileType = sourcecode.cpp.cpp; name = Viscosity.cpp; path = Parameter/Solid/Viscosity.cpp; sourceTree = "<group>"; };
BEE934761CFD8B4F00158440 /* MatrixConversion.cpp */ = {isa = PBXFileReference; fileEncoding = 4; lastKnownFileType = sourcecode.cpp.cpp; name = MatrixConversion.cpp; path = ConvertLinearAlgebra/Freefem/MatrixConversion.cpp; sourceTree = "<group>"; };
BEE934771CFD8B4F00158440 /* MatrixConversion.hpp */ = {isa = PBXFileReference; fileEncoding = 4; lastKnownFileType = sourcecode.cpp.h; name = MatrixConversion.hpp; path = ConvertLinearAlgebra/Freefem/MatrixConversion.hpp; sourceTree = "<group>"; };
......@@ -8507,6 +8510,7 @@
BE7E68842069554D00AA2FB3 /* TupleHasType */ = {
isa = PBXGroup;
children = (
BEE79C6C2136F51E00388FED /* CMakeLists.txt */,
BE7E6887206955AD00AA2FB3 /* main.cpp */,
);
path = TupleHasType;
......@@ -10122,6 +10126,7 @@
children = (
BE8B5FA62077D65F00DC005E /* CMakeLists.txt */,
BE8B5F8D207792B600DC005E /* Tclap */,
BEE79C672136F43000388FED /* NanOrInf */,
BEC5CFB41F4DBF9B00A2863B /* Mpi */,
);
path = ThirdParty;
......@@ -10410,6 +10415,15 @@
name = Products;
sourceTree = "<group>";
};
BEE79C672136F43000388FED /* NanOrInf */ = {
isa = PBXGroup;
children = (
BEE79C6A2136F43000388FED /* CMakeLists.txt */,
BEE79C6B2136F43000388FED /* test.cpp */,
);
path = NanOrInf;
sourceTree = "<group>";
};
BEE9346E1CFD8ACB00158440 /* ConvertLinearAlgebra */ = {
isa = PBXGroup;
children = (
Sources/FiniteElement/FiniteElementSpace/GodOfDof.cpp
View file @ 58de6df4
......@@ -1159,20 +1159,7 @@ namespace MoReFEM
return;
const auto& value_list = dof_storage.GetDofValueList();
# ifdef MOREFEM_CHECK_NAN_AND_INF
if (std::any_of(value_list.cbegin(),
value_list.cend(),
[](double val)
{
return std::isnan(val) || std::isinf(val);
}))
{
throw Exception("Nan or Inf value found in GodOfDof::ApplyPseudoElimination().", __FILE__, __LINE__);
}
# endif // MOREFEM_CHECK_NAN_AND_INF
vector.SetValues(dof_index_list,
value_list.data(),
INSERT_VALUES,
......
Sources/FormulationSolver/VariationalFormulation.hxx
View file @ 58de6df4
......@@ -112,10 +112,6 @@ namespace MoReFEM
{
assert(out.GetNumberingSubset() == rhs.GetNumberingSubset());
#ifdef MOREFEM_CHECK_NAN_AND_INF
Wrappers::Petsc::CheckNumericValues(this->GetMpi(), "rhs", rhs, invoking_file, invoking_line);
#endif // MOREFEM_CHECK_NAN_AND_INF
auto do_print_solver_info =
(this->GetTimeManager().NtimeModified() % GetDisplayValue()) == 0
? Wrappers::Petsc::print_solver_infos::yes
......@@ -153,11 +149,6 @@ namespace MoReFEM
}
out.UpdateGhosts(__FILE__, __LINE__);
#ifdef MOREFEM_CHECK_NAN_AND_INF
Wrappers::Petsc::CheckNumericValues(this->GetMpi(), "out", out, invoking_file, invoking_line);
#endif // MOREFEM_CHECK_NAN_AND_INF
}
......
Sources/OperatorInstances/VariationalOperator/BilinearForm/Local/ScalarDivVectorial.cpp
View file @ 58de6df4
......@@ -11,10 +11,6 @@
#include "ThirdParty/Wrappers/Seldon/SeldonFunctions.hpp"
#ifdef MOREFEM_CHECK_NAN_AND_INF
# include "ThirdParty/Wrappers/Seldon/MatrixOperations.hpp"
#endif // MOREFEM_CHECK_NAN_AND_INF
#include "OperatorInstances/VariationalOperator/BilinearForm/Local/ScalarDivVectorial.hpp"
......
Sources/OperatorInstances/VariationalOperator/NonlinearForm/Local/SecondPiolaKirchhoffStressTensor.hpp
View file @ 58de6df4
......@@ -20,12 +20,6 @@
# include "Utilities/LinearAlgebra/Storage/Local/LocalMatrixStorage.hpp"
# include "Utilities/LinearAlgebra/Storage/Local/LocalVectorStorage.hpp"
# ifdef MOREFEM_CHECK_NAN_AND_INF
# include "ThirdParty/Wrappers/Seldon/MatrixOperations.hpp"
# endif // MOREFEM_CHECK_NAN_AND_INF
# include "Parameters/InitParameterFromInputData/InitParameterFromInputData.hpp"
# include "ParameterInstances/Compound/Solid/Solid.hpp"
......
Sources/OperatorInstances/VariationalOperator/NonlinearForm/Local/SecondPiolaKirchhoffStressTensor.hxx
View file @ 58de6df4
......@@ -476,10 +476,6 @@ namespace MoReFEM
dW,
d2W);
# ifdef MOREFEM_CHECK_NAN_AND_INF
Wrappers::Seldon::ThrowIfNanInside(d2W, __FILE__, __LINE__);
# endif // MOREFEM_CHECK_NAN_AND_INF
// ===================================================================================
// Finally build the terms that are actually required.
// ===================================================================================
......
Sources/OperatorInstances/VariationalOperator/NonlinearForm/Local/SecondPiolaKirchhoffStressTensor/HyperelasticityPolicy/Hyperelasticity.hxx
View file @ 58de6df4
......@@ -87,6 +87,10 @@ namespace MoReFEM
const auto& dI3dC =
invariant_holder.GetFirstDerivativeWrtCauchyGreen(invariant_holder_type::invariants_first_derivative_index::dI3dC);
// Note: one of these quantities might be infinite... but in this case this infinite (or a resulting
// nan) will be assembled into a Petsc matrix or vector, which will in turn handle properly
// the case. See #1317 to know more about it.
const double dWdI1 = law.FirstDerivativeWFirstInvariant(invariant_holder, quad_pt, geom_elt);
const double dWdI2 = law.FirstDerivativeWSecondInvariant(invariant_holder, quad_pt, geom_elt);
const double dWdI3 = law.FirstDerivativeWThirdInvariant(invariant_holder, quad_pt, geom_elt);
......@@ -100,12 +104,6 @@ namespace MoReFEM
// d2I1dCdC = d2I4dCdC = 0
Seldon::Add(dWdI2, d2I2dCdC, d2W);
Seldon::Add(dWdI3, d2I3dCdC, d2W);
# ifdef MOREFEM_CHECK_NAN_AND_INF
Wrappers::Seldon::ThrowIfNanInside(d2I2dCdC, __FILE__, __LINE__);
Wrappers::Seldon::ThrowIfNanInside(d2I3dCdC, __FILE__, __LINE__);
Wrappers::Seldon::ThrowIfNanInside(d2W, __FILE__, __LINE__);
# endif // MOREFEM_CHECK_NAN_AND_INF
}
{
......
Sources/OperatorInstances/VariationalOperator/NonlinearForm/Local/SecondPiolaKirchhoffStressTensor/Internal/Helper.cpp
View file @ 58de6df4
......@@ -9,10 +9,6 @@
/// \addtogroup OperatorInstancesGroup
/// \{
#ifdef MOREFEM_CHECK_NAN_AND_INF
# include "ThirdParty/Wrappers/Seldon/MatrixOperations.hpp"
#endif // MOREFEM_CHECK_NAN_AND_INF
#include "OperatorInstances/VariationalOperator/NonlinearForm/Local/SecondPiolaKirchhoffStressTensor/Internal/Helper.hpp"
......@@ -81,10 +77,6 @@ namespace MoReFEM
{
Seldon::Mlt(1., transposed_De, d2W, linear_part_intermediate_matrix);
Seldon::Mlt(1., linear_part_intermediate_matrix, De, linear_part);
# ifdef MOREFEM_CHECK_NAN_AND_INF
::MoReFEM::Wrappers::Seldon::ThrowIfNanInside(linear_part, __FILE__, __LINE__);
# endif // MOREFEM_CHECK_NAN_AND_INF
}
......
Sources/Operators/ConformInterpolator/Internal/ComputePatternHelper.cpp
View file @ 58de6df4
......@@ -123,10 +123,6 @@ namespace MoReFEM
const auto& source_unknown_storage = source_data.GetExtendedUnknownList();
const auto& target_unknown_storage = target_data.GetExtendedUnknownList();
# ifdef MOREFEM_CHECK_NAN_AND_INF
assert(!::MoReFEM::Wrappers::Seldon::IsZeroMatrix(local_projection_matrix));
# endif // MOREFEM_CHECK_NAN_AND_INF
const auto interpreted_local_projection_matrix = LocalMatrixContent(local_projection_matrix);
const auto Nlocal_row = static_cast<std::size_t>(local_projection_matrix.GetM());
......
Sources/Test/ThirdParty/CMakeLists.txt
View file @ 58de6df4
include(${CMAKE_CURRENT_LIST_DIR}/Mpi/CMakeLists.txt)
include(${CMAKE_CURRENT_LIST_DIR}/Tclap/CMakeLists.txt)
include(${CMAKE_CURRENT_LIST_DIR}/NanOrInf/CMakeLists.txt)
\ No newline at end of file
Sources/Test/ThirdParty/NanOrInf/CMakeLists.txt 0 → 100644
View file @ 58de6df4
add_executable(MoReFEMTestNanOrInf
${CMAKE_CURRENT_LIST_DIR}/test.cpp
)
target_link_libraries(MoReFEMTestNanOrInf
${MOREFEM_UTILITIES})
add_test(NanOrInf MoReFEMTestNanOrInf)
morefem_install(MoReFEMTestNanOrInf)
Sources/Test/ThirdParty/NanOrInf/test.cpp 0 → 100644
View file @ 58de6df4
//! \file
//
//
#define CATCH_CONFIG_MAIN
#include "ThirdParty/Source/Catch/catch.hpp"
#include "ThirdParty/Wrappers/Seldon/MatrixOperations.hpp"
using namespace MoReFEM;
namespace // anonymous
{
double local_inf = 1. / 0.;
double local_nan = local_inf * 0.;
}
TEST_CASE("Check IsNumber() is behaving as expected.")
{
CHECK(NumericNS::IsNumber(5.));
CHECK(NumericNS::IsNumber(4.f));
CHECK(NumericNS::IsNumber(local_inf) == false);
CHECK(NumericNS::IsNumber(local_nan) == false);
}
TEST_CASE("Check the function to check a Seldon vector.")
{
LocalVector vector(3);
vector(0) = 0.;
vector(1) = -10.;
vector(2) = 10.;
CHECK(Wrappers::Seldon::AreNanOrInfValues(vector) == false);
vector(2) = local_nan;
CHECK(Wrappers::Seldon::AreNanOrInfValues(vector) == true);
vector(2) = 0.;
vector(0) = local_inf;
CHECK(Wrappers::Seldon::AreNanOrInfValues(vector) == true);
}
TEST_CASE("Check the function to check a Seldon matrix.")
{
LocalMatrix matrix(3, 3);
matrix.Zero();
CHECK(Wrappers::Seldon::AreNanOrInfValues(matrix) == false);
matrix(2, 0) = local_nan;
CHECK(Wrappers::Seldon::AreNanOrInfValues(matrix) == true);
matrix(2, 0) = 0.;
matrix(1, 1) = local_inf;
CHECK(Wrappers::Seldon::AreNanOrInfValues(matrix) == true);
}
Sources/ThirdParty/Wrappers/Petsc/Vector/Vector.cpp
View file @ 58de6df4
......@@ -681,42 +681,7 @@ namespace MoReFEM
throw ExceptionNS::Exception(error_code, "VecRestoreArrayRead", invoking_file, invoking_line);
}
#ifdef MOREFEM_CHECK_NAN_AND_INF
void CheckNumericValues(const Mpi& mpi,
const std::string& vector_name,
const Vector& vector,
const char* invoking_file, int invoking_line)
{
AccessVectorContent<Utilities::Access::read_only> content(vector, invoking_file, invoking_line);
const auto size = content.GetSize(invoking_file, invoking_line);
for (unsigned int i = 0u; i < size; ++i)
{
const auto current_value = content.GetValue(i);
if (std::isnan(current_value))
{
std::ostringstream oconv;
oconv << "Nan value found in vector " << vector_name << " called on " << invoking_file << ", line "
<< invoking_line << '.';
::MoReFEM::ExceptionNS::PrintAndAbort(mpi, oconv.str());
}
if (std::isinf(current_value))
{
std::ostringstream oconv;
oconv << "Inf value found in vector " << vector_name << " called on " << invoking_file << ", line "
<< invoking_line << '.';
::MoReFEM::ExceptionNS::PrintAndAbort(mpi, oconv.str());
}
}
}
#endif // MOREFEM_CHECK_NAN_AND_INF
#ifdef MOREFEM_CHECK_UPDATE_GHOSTS_CALL_RELEVANCE
......
Sources/ThirdParty/Wrappers/Petsc/Vector/Vector.hpp
View file @ 58de6df4
......@@ -820,26 +820,6 @@ namespace MoReFEM
// ============================
# ifdef MOREFEM_CHECK_NAN_AND_INF
/*!
* \brief Check none of the values is inf or nan.
*
* \copydetails doxygen_hide_mpi_param
* \param[in] vector_name Name of the vector under check; it is only used to provide information in the log.
* \param[in] vector Vector being scrutinized.
* \copydoc doxygen_hide_invoking_file_and_line
*/
void CheckNumericValues(const Mpi& mpi,
const std::string& vector_name,
const Vector& vector,
const char* invoking_file, int invoking_line);
# endif // MOREFEM_CHECK_NAN_AND_INF
} // namespace Petsc
......
Sources/ThirdParty/Wrappers/Seldon/MatrixOperations.cpp
View file @ 58de6df4
......@@ -9,6 +9,8 @@
/// \addtogroup ThirdPartyGroup
/// \{
#include <array>
#include "ThirdParty/Wrappers/Seldon/MatrixOperations.hpp"
#include "Utilities/Numeric/Numeric.hpp"
......@@ -97,33 +99,7 @@ namespace MoReFEM
}
}
#ifdef MOREFEM_CHECK_NAN_AND_INF
void ThrowIfNanInside(const LocalMatrix& matrix, const char* invoking_file, int invoking_line)
{
const int Nrow = matrix.GetM();
const int Ncol = matrix.GetN();
assert(Nrow > 0);
assert(Ncol > 0);
bool is_nan = false;
for (int i = 0; !is_nan && i < Nrow; ++i)
{
for (int j = 0; !is_nan && j < Ncol; ++j)
{
if (std::isnan(matrix(i,j)))
is_nan = true;
}
}
if (is_nan)
throw Exception("NaN value found in local matrix!", invoking_file, invoking_line);
}
#endif // MOREFEM_CHECK_NAN_AND_INF
bool IsZeroMatrix(const LocalMatrix& matrix)
{
const auto Nrow = matrix.GetM();
......
Sources/ThirdParty/Wrappers/Seldon/MatrixOperations.hpp
View file @ 58de6df4
......@@ -107,30 +107,38 @@ namespace MoReFEM
#endif // NDEBUG
#ifdef MOREFEM_CHECK_NAN_AND_INF
/*!
* \brief Check whether a matrix is filled only with zero.
*/
bool IsZeroMatrix(const LocalMatrix& matrix);
/*!
* \brief Throw an exception if there is a nan values in a matrix.
* \brief Check whether there is a nan or inf value in a Seldon matrix.
*
* \param[in] matrix Matrix being checked.
* \param[in] invoking_file Usually the macro __FILE__; the point is to identify where the exception was
* thrown first.
* \param[in] invoking_line Usually the macro __LINE__.
* \tparam MatrixT A dense Seldon matrix.
*
* This method is first and foremost for debug mode, as there is an obvious efficiency cost (the function
* is O(n), as std::isnan is checked upon every value of the matrix).
* \param[in] matrix Matrix being scrutinized.
*
* \return True if at least one value if NaN or inf.
*/
void ThrowIfNanInside(const LocalMatrix& matrix, const char* invoking_file, int invoking_line);
#endif // MOREFEM_CHECK_NAN_AND_INF
template<class MatrixT>
std::enable_if_t<Utilities::IsSpecializationOf<::Seldon::Matrix, MatrixT>::value, bool>
AreNanOrInfValues(const MatrixT& matrix);
/*!
* \brief Check whether a matrix is filled only with zero.
* \brief Check whether there is a nan or inf value in a Seldon vector.
*
* \tparam VectorT A dense Seldon vector.
*
* \param[in] vector Vector being scrutinized.
*
* \return True if at least one value if NaN or inf.
*/
bool IsZeroMatrix(const LocalMatrix& matrix);
template<class VectorT>
std::enable_if_t<Utilities::IsSpecializationOf<::Seldon::Vector, VectorT>::value, bool>
AreNanOrInfValues(const VectorT& vector);
} // namespace Seldon
......
Sources/ThirdParty/Wrappers/Seldon/MatrixOperations.hxx
View file @ 58de6df4
......@@ -79,6 +79,46 @@ namespace MoReFEM
#endif // NDEBUG
template<class MatrixT>
std::enable_if_t<Utilities::IsSpecializationOf<::Seldon::Matrix, MatrixT>::value, bool>
AreNanOrInfValues(const MatrixT& matrix)
{
const int Nrow = matrix.GetM();
const int Ncol = matrix.GetN();
assert(Nrow > 0);
assert(Ncol > 0);
for (int i = 0; i < Nrow; ++i)
{
for (int j = 0; j < Ncol; ++j)
{
if (!NumericNS::IsNumber(matrix(i, j)))
return true;
}
}
return false;
}
template<class VectorT>
std::enable_if_t<Utilities::IsSpecializationOf<::Seldon::Vector, VectorT>::value, bool>
AreNanOrInfValues(const VectorT& vector)
{
const int size = vector.GetSize();
assert(size > 0);
for (int i = 0; i < size; ++i)
{
if (!NumericNS::IsNumber(vector(i)))
return true;
}
return false;
}
} // namespace Seldon
......
Sources/Utilities/Numeric/Numeric.hpp
View file @ 58de6df4
......@@ -147,6 +147,18 @@ namespace MoReFEM
AreEqual(T lhs, T rhs, T epsilon = DefaultEpsilon<T>()) noexcept;
/*!
* \brief Checks whether a quantity is finite and defined, and if not say so explicitly.
*
* \param[in] value The value being addressed.
*
* \return True if the value is neither inf nor nan.
*/
template<class T>
std::enable_if_t<std::is_floating_point<T>::value, bool>
IsNumber(T value) noexcept;
} // namespace NumericNS
......
Sources/Utilities/Numeric/Numeric.hxx
View file @ 58de6df4
......@@ -134,6 +134,14 @@ namespace MoReFEM
}
template<class T>
std::enable_if_t<std::is_floating_point<T>::value, bool> IsNumber(T value) noexcept
{
return (std::isnan(value) || std::isinf(value)) ? false : true;
}
} // namespace NumericNS
......
Markdown is supported
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment