Commit e7c879ab authored by GILLES Sebastien's avatar GILLES Sebastien
Browse files

#890 Core: complete implementation of AssertMatrixRespectPattern() and test it in poromechanics.

parent 79637520
......@@ -8,6 +8,13 @@
// Copyright (c) 2015 Inria. All rights reserved.
//
#include <algorithm>
#include "Utilities/Containers/Print.hpp" //\todo #820 DEV
#include "Utilities/Exceptions/Exception.hpp"
#include "Utilities/Exceptions/PrintAndAbort.hpp"
#include "ThirdParty/Wrappers/Petsc/Matrix/MatrixPattern.hpp"
#include "Core/LinearAlgebra/GlobalMatrix.hpp"
#include "Core/NumberingSubset.hpp"
......@@ -41,5 +48,69 @@ namespace HappyHeart
}
#ifndef NDEBUG
void AssertMatrixRespectPattern(const Wrappers::Mpi& mpi,
const GlobalMatrix& matrix,
const Wrappers::Petsc::MatrixPattern& pattern,
const char* invoking_file, int invoking_line)
{
try
{
const auto Nrow = matrix.GetProcessorWiseSize(__FILE__, __LINE__).first;
if (Nrow != static_cast<int>(pattern.Nrow()))
throw Exception("Number of rows in pattern and matrix is not the same!",
invoking_file, invoking_line);
decltype(auto) iCsr = pattern.GetICsr();
decltype(auto) jCsr = pattern.GetJCsr();
assert(iCsr.size() == static_cast<std::size_t>(Nrow + 1)
&& "This would highlight a bug within MatrixPattern class");
Utilities::PrintContainer(jCsr, std::cout, ", ", "Assert jCSR -> [");
auto current_jcsr_index = 0ul;
for (auto i = 0; i < Nrow; ++i)
{
std::vector<PetscInt> position_list_in_matrix;
std::vector<PetscScalar> value_list_in_matrix;
matrix.GetRow(i, position_list_in_matrix, value_list_in_matrix, __FILE__, __LINE__);
const auto Nvalue_in_pattern_row =
static_cast<std::size_t>(iCsr[static_cast<std::size_t>(i + 1)] - iCsr[static_cast<std::size_t>(i)]);
std::vector<PetscInt> position_list_in_pattern;
for (auto j = 0ul; j < Nvalue_in_pattern_row; ++j)
{
assert(current_jcsr_index < jCsr.size());
position_list_in_pattern.push_back(static_cast<PetscInt>(jCsr[current_jcsr_index + j]));
}
current_jcsr_index += Nvalue_in_pattern_row;
Utilities::PrintContainer(position_list_in_pattern, std::cout, ", ", "Pattern -> [");
Utilities::PrintContainer(position_list_in_matrix, std::cout, ", ", "Matrix -> [");
if (!std::includes(position_list_in_pattern.cbegin(),
position_list_in_pattern.cend(),
position_list_in_matrix.cbegin(),
position_list_in_matrix.cend()))
{
std::cout << "EXCEPTION!!!!" << std::endl;
throw Exception("Position in the actual matrix must match those defined in the pattern!",
invoking_file, invoking_line);
}
}
}
catch(const Exception& e)
{
ExceptionNS::PrintAndAbort(mpi, e.what());
}
}
#endif // NDEBUG
} // namespace HappyHeart
......@@ -97,6 +97,25 @@ namespace HappyHeart
};
# ifndef NDEBUG
/*!
* \brief Check whether a matrix really follows the pattern expected from its numbering subset.
*
* If not fulfilled, the program is aborted.
*
* \copydetails doxygen_hide_mpi_param
* \param[in] matrix Matrix being investigated.
* \param[in] pattern Object that describe the sparse pattern the matrix is expected to follow. In HappyHeart
* such objects are stored in the GodOfDof.
* \copydetails doxygen_hide_invoking_file_and_line
*/
void AssertMatrixRespectPattern(const Wrappers::Mpi& mpi,
const GlobalMatrix& matrix,
const Wrappers::Petsc::MatrixPattern& pattern,
const char* invoking_file, int invoking_line);
# endif // NDEBUG
/*!
* \brief Swap two matrices.
......
......@@ -10,6 +10,8 @@
#include "ModelInstances/UnderDevelopment/Poromechanics/ImplicitStepFluidVariationalFormulation.hpp"
#include "FiniteElement/FiniteElementSpace/GodOfDofManager.hpp" // \todo #820 For DEV only!
namespace HappyHeart
{
......@@ -25,6 +27,8 @@ namespace HappyHeart
void VariationalFormulation::Perform()
{
decltype(auto) mpi = parent::MpiHappyHeart();
decltype(auto) time_manager = parent::GetTimeManager();
const auto time = time_manager.GetTime();
decltype(auto) god_of_dof = parent::GetGodOfDof();
......@@ -76,6 +80,35 @@ namespace HappyHeart
GlobalMatrixWithCoefficient fluid_matrix_on_solid_mesh(fluid_mass_matrix_on_solid_mesh, 1. / time_step);
GetMassOperatorFluidMassOnSolidMesh().Assemble(std::make_tuple(std::ref(fluid_matrix_on_solid_mesh)));
fluid_mass_matrix_on_solid_mesh.View(MpiHappyHeart(), __FILE__, __LINE__);
{
decltype(auto) solid_god_of_dof = GodOfDofManager::GetInstance().GetGodOfDof(2);
decltype(auto) pattern = solid_god_of_dof.GetMatrixPattern(fluid_mass_matrix_on_solid_mesh.GetRowNumberingSubset(),
fluid_mass_matrix_on_solid_mesh.GetRowNumberingSubset());
decltype(auto) icsr = pattern.GetICsr();
decltype(auto) jcsr = pattern.GetJCsr();
{
std::ostringstream oconv;
oconv << MpiHappyHeart().GetRankPreffix();
Utilities::PrintContainer(icsr, oconv, ", ", "iCSR = [");
Utilities::PrintContainer(jcsr, oconv, ", ", "jCSR = [");
std::cout << oconv.str();
}
#ifndef NDEBUG
AssertMatrixRespectPattern(mpi,
fluid_mass_matrix_on_solid_mesh,
pattern,
__FILE__, __LINE__);
#endif // NDEBUG
}
{
auto& t11 = GetNonCstFluidMassMatrixOnFluidMesh();
......@@ -84,6 +117,18 @@ namespace HappyHeart
t11,
__FILE__, __LINE__);
#ifndef NDEBUG
{
decltype(auto) fluid_mass_numbering_subset =
god_of_dof.GetNumberingSubset(EnumUnderlyingType(NumberingSubsetIndex::fluid_mass));
decltype(auto) pattern = god_of_dof.GetMatrixPattern(fluid_mass_numbering_subset,
fluid_mass_numbering_subset);
AssertMatrixRespectPattern(mpi, t11, pattern, __FILE__, __LINE__);
}
#endif // NDEBUG
t11.View(parent::MpiHappyHeart(),
parent::GetOutputDirectory(numbering_subset) + "/T11_" + std::to_string(time_manager.NtimeModified()) + ".m",
__FILE__, __LINE__,
......
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