Commit 9fd2c1bd authored by GILLES Sebastien's avatar GILLES Sebastien
Browse files

#859 Petsc: add MatTransposeMatMult Symbolic/Numeric (through a hack as they...

#859 Petsc: add MatTransposeMatMult Symbolic/Numeric (through a hack as they don't exist in Petsc...)
parent e7480106
...@@ -21,7 +21,7 @@ transient = { ...@@ -21,7 +21,7 @@ transient = {
-- Maximum time, if set to zero run a static case. -- Maximum time, if set to zero run a static case.
-- Expected format: VALUE -- Expected format: VALUE
-- Constraint: v >= 0. -- Constraint: v >= 0.
timeMax = 2.e-3 timeMax = 1.e-3
} -- transient } -- transient
......
...@@ -80,53 +80,66 @@ namespace HappyHeart ...@@ -80,53 +80,66 @@ namespace HappyHeart
GlobalMatrixWithCoefficient fluid_matrix_on_solid_mesh(fluid_mass_matrix_on_solid_mesh, 1. / time_step); 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))); GetMassOperatorFluidMassOnSolidMesh().Assemble(std::make_tuple(std::ref(fluid_matrix_on_solid_mesh)));
decltype(auto) solid_god_of_dof = GodOfDofManager::GetInstance().GetGodOfDof(2);
AssertMatrixRespectPattern(solid_god_of_dof,
fluid_mass_matrix_on_solid_mesh,
__FILE__, __LINE__);
fluid_mass_matrix_on_solid_mesh.View(MpiHappyHeart(), __FILE__, __LINE__); fluid_mass_matrix_on_solid_mesh.View(MpiHappyHeart(), __FILE__, __LINE__);
{ {
decltype(auto) solid_god_of_dof = GodOfDofManager::GetInstance().GetGodOfDof(2); auto& t11 = GetNonCstFluidMassMatrixOnFluidMesh();
// \todo #820 Clearly first approach here:
// - If with help of Petsc dev I manage to make MatPtAP, use that instead!
// - Symbolic/Numeric should be used as well...
// Currently we decompose computation of tPAP:
auto& intermediate = GetNonCstWorkMatrixFluidMassMatrixOnFluidMesh();
decltype(auto) pattern = solid_god_of_dof.GetMatrixPattern(fluid_mass_matrix_on_solid_mesh.GetRowNumberingSubset(), Wrappers::Petsc::MatMatMult(fluid_mass_matrix_on_solid_mesh,
fluid_mass_matrix_on_solid_mesh.GetRowNumberingSubset()); inverse_interp_mass_from_solid_mesh_->GetInterpolationMatrix(),
intermediate,
__FILE__, __LINE__);
Wrappers::Petsc::MatTransposeMatMult(inverse_interp_mass_from_solid_mesh_->GetInterpolationMatrix(),
intermediate,
t11,
__FILE__, __LINE__);
// Wrappers::Petsc::PtAP(fluid_mass_matrix_on_solid_mesh,
// GetInterpMassFromSolidMesh().GetInterpolationMatrix(),
// foo,
// __FILE__, __LINE__);
// foo.View(MpiHappyHeart(), __FILE__, __LINE__);
{
decltype(auto) pattern = god_of_dof.GetMatrixPattern(t11.GetRowNumberingSubset(),
t11.GetColNumberingSubset());
decltype(auto) icsr = pattern.GetICsr(); decltype(auto) icsr = pattern.GetICsr();
decltype(auto) jcsr = pattern.GetJCsr(); decltype(auto) jcsr = pattern.GetJCsr();
{ {
std::ostringstream oconv; std::ostringstream oconv;
oconv << MpiHappyHeart().GetRankPreffix(); oconv << MpiHappyHeart().GetRankPreffix();
Utilities::PrintContainer(icsr, oconv, ", ", "iCSR = ["); Utilities::PrintContainer(icsr, oconv, ", ", "iCSR = [");
Utilities::PrintContainer(jcsr, oconv, ", ", "jCSR = ["); Utilities::PrintContainer(jcsr, oconv, ", ", "jCSR = [");
std::cout << oconv.str(); std::cout << oconv.str();
} }
#ifndef NDEBUG #ifndef NDEBUG
AssertMatrixRespectPattern(mpi, AssertMatrixRespectPattern(god_of_dof,
fluid_mass_matrix_on_solid_mesh, t11,
pattern,
__FILE__, __LINE__); __FILE__, __LINE__);
#endif // NDEBUG #endif // NDEBUG
} }
{
auto& t11 = GetNonCstFluidMassMatrixOnFluidMesh();
Wrappers::Petsc::MatMatMultNumeric(GetInterpMassFromSolidMesh().GetInterpolationMatrix(),
fluid_mass_matrix_on_solid_mesh,
t11,
__FILE__, __LINE__);
#ifndef NDEBUG #ifndef NDEBUG
{ AssertMatrixRespectPattern(god_of_dof, t11, __FILE__, __LINE__);
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 #endif // NDEBUG
t11.View(parent::MpiHappyHeart(), t11.View(parent::MpiHappyHeart(),
......
...@@ -314,6 +314,12 @@ namespace HappyHeart ...@@ -314,6 +314,12 @@ namespace HappyHeart
//! Non constant accessor to the work matrix that gets the exact same structure as system matrix. //! Non constant accessor to the work matrix that gets the exact same structure as system matrix.
GlobalMatrix& GetNonCstWorkMatrixSystemLike() noexcept; GlobalMatrix& GetNonCstWorkMatrixSystemLike() noexcept;
/*!
* \brief Non constant accessor to the work matrix that share the pattern of
* fluid_mass_matrix_on_fluid_mesh_. There only as long as MatPtAP doesn't work.
*/
GlobalMatrix& GetNonCstWorkMatrixFluidMassMatrixOnFluidMesh() noexcept;
///@} ///@}
...@@ -344,6 +350,8 @@ namespace HappyHeart ...@@ -344,6 +350,8 @@ namespace HappyHeart
//! Interpolator to obtain on fluid dofs the result of mass operator applied on solid mesh, //! Interpolator to obtain on fluid dofs the result of mass operator applied on solid mesh,
NonConformInterpolatorNS::FromVertexMatching::const_unique_ptr interp_mass_from_solid_mesh_ = nullptr; NonConformInterpolatorNS::FromVertexMatching::const_unique_ptr interp_mass_from_solid_mesh_ = nullptr;
NonConformInterpolatorNS::FromVertexMatching::const_unique_ptr inverse_interp_mass_from_solid_mesh_ = nullptr;
//! Interpolator to expand from fluid mass to monolithic. //! Interpolator to expand from fluid mass to monolithic.
ConformInterpolatorNS::SubsetOrSuperset::unique_ptr expand_mass_ = nullptr; ConformInterpolatorNS::SubsetOrSuperset::unique_ptr expand_mass_ = nullptr;
...@@ -390,7 +398,6 @@ namespace HappyHeart ...@@ -390,7 +398,6 @@ namespace HappyHeart
GlobalMatrix::unique_ptr fluid_mass_matrix_on_solid_mesh_ = nullptr; GlobalMatrix::unique_ptr fluid_mass_matrix_on_solid_mesh_ = nullptr;
//! Matrix obtained from \a fluid_mass_matrix_on_solid_mesh_ following interpolation. //! Matrix obtained from \a fluid_mass_matrix_on_solid_mesh_ following interpolation.
GlobalMatrix::unique_ptr fluid_mass_matrix_on_fluid_mesh_ = nullptr; GlobalMatrix::unique_ptr fluid_mass_matrix_on_fluid_mesh_ = nullptr;
//! Interpolation matrix to go from fluid mass on solid mesh to the monolithic matrix on fluid mesh. //! Interpolation matrix to go from fluid mass on solid mesh to the monolithic matrix on fluid mesh.
...@@ -402,6 +409,12 @@ namespace HappyHeart ...@@ -402,6 +409,12 @@ namespace HappyHeart
//! Work matrix that gets the exact same structure as system matrix. //! Work matrix that gets the exact same structure as system matrix.
GlobalMatrix::unique_ptr work_matrix_system_like_ = nullptr; GlobalMatrix::unique_ptr work_matrix_system_like_ = nullptr;
/*!
* \brief Work matrix that share the pattern of fluid_mass_matrix_on_fluid_mesh_. There only as long as
* MatPtAP doesn't work.
*/
GlobalMatrix::unique_ptr work_matrix_fluid_mass_matrix_on_fluid_mesh_ = nullptr;
///@} ///@}
......
...@@ -291,7 +291,13 @@ namespace HappyHeart ...@@ -291,7 +291,13 @@ namespace HappyHeart
{ {
assert(!(!work_matrix_system_like_)); assert(!(!work_matrix_system_like_));
return *work_matrix_system_like_; return *work_matrix_system_like_;
}
inline GlobalMatrix& VariationalFormulation::GetNonCstWorkMatrixFluidMassMatrixOnFluidMesh() noexcept
{
assert(!(!work_matrix_fluid_mass_matrix_on_fluid_mesh_));
return *work_matrix_fluid_mass_matrix_on_fluid_mesh_;
} }
......
...@@ -106,6 +106,8 @@ namespace HappyHeart ...@@ -106,6 +106,8 @@ namespace HappyHeart
std::make_unique<GlobalMatrix>(fluid_mass_numbering_subset, fluid_mass_numbering_subset); std::make_unique<GlobalMatrix>(fluid_mass_numbering_subset, fluid_mass_numbering_subset);
work_matrix_fluid_mass_matrix_on_fluid_mesh_ =
std::make_unique<GlobalMatrix>(fluid_mass_numbering_subset, fluid_mass_numbering_subset);
} }
...@@ -211,6 +213,13 @@ namespace HappyHeart ...@@ -211,6 +213,13 @@ namespace HappyHeart
std::make_unique<type>(input_parameter_data, std::make_unique<type>(input_parameter_data,
EnumUnderlyingType(InitVertexMatchingInterpolator::solid_mass), EnumUnderlyingType(InitVertexMatchingInterpolator::solid_mass),
EnumUnderlyingType(InitVertexMatchingInterpolator::fluid_mass)); EnumUnderlyingType(InitVertexMatchingInterpolator::fluid_mass));
inverse_interp_mass_from_solid_mesh_ =
std::make_unique<type>(input_parameter_data,
EnumUnderlyingType(InitVertexMatchingInterpolator::fluid_mass),
EnumUnderlyingType(InitVertexMatchingInterpolator::solid_mass));
} }
{ {
...@@ -231,7 +240,6 @@ namespace HappyHeart ...@@ -231,7 +240,6 @@ namespace HappyHeart
mesh_dimension, mesh_dimension,
default_quadrature_rule_set, default_quadrature_rule_set,
DoComputeProcessorWiseLocal2Global::yes); // \todo #820 Check if no is ok. DoComputeProcessorWiseLocal2Global::yes); // \todo #820 Check if no is ok.
} }
...@@ -254,7 +262,6 @@ namespace HappyHeart ...@@ -254,7 +262,6 @@ namespace HappyHeart
god_of_dof.GetNumberingSubset(EnumUnderlyingType(NumberingSubsetIndex::fluid_mass))); god_of_dof.GetNumberingSubset(EnumUnderlyingType(NumberingSubsetIndex::fluid_mass)));
reduce_to_mass_->Init(); reduce_to_mass_->Init();
} }
...@@ -278,36 +285,37 @@ namespace HappyHeart ...@@ -278,36 +285,37 @@ namespace HappyHeart
} }
{
auto& interpolator = GetNonCstFluidMassInterpolationMatrix();
const auto& foo = GetInterpMassFromSolidMesh().GetInterpolationMatrix();
Wrappers::Petsc::MatMatMult(GetExpandMass().GetInterpolationMatrix(),
foo,
interpolator,
__FILE__, __LINE__);
auto& inverse_interpolator = GetNonCstInverseFluidMassInterpolationMatrix();
GlobalMatrix transposed(foo.GetColNumberingSubset(), foo.GetRowNumberingSubset());
Wrappers::Petsc::MatTransposeMatMult(foo,
GetReduceToMass().GetInterpolationMatrix(),
inverse_interpolator,
__FILE__, __LINE__);
}
{ {
// Allocate the pattern of a matrix from the result of the product of two other matrices. // Allocate the pattern of a matrix from the result of the product of two other matrices.
auto& fluid_mass_matrix_on_fluid_mesh = GetNonCstFluidMassMatrixOnFluidMesh();
const auto& fluid_mass_matrix_on_solid_mesh = GetFluidMassMatrixOnSolidMesh(); const auto& fluid_mass_matrix_on_solid_mesh = GetFluidMassMatrixOnSolidMesh();
Wrappers::Petsc::MatMatMultSymbolic(GetInterpMassFromSolidMesh().GetInterpolationMatrix(), auto& t11 = GetNonCstFluidMassMatrixOnFluidMesh();
fluid_mass_matrix_on_solid_mesh, auto& intermediate = GetNonCstWorkMatrixFluidMassMatrixOnFluidMesh();
fluid_mass_matrix_on_fluid_mesh,
__FILE__, __LINE__); // decltype(auto) interp_matrix = GetInterpMassFromSolidMesh().GetInterpolationMatrix();
//
// decltype(auto) size1 = fluid_mass_matrix_on_solid_mesh.GetProcessorWiseSize(__FILE__, __LINE__);
// decltype(auto) size2 = interp_matrix.GetProcessorWiseSize(__FILE__, __LINE__);
// // decltype(auto) size3 = intermediate.GetProcessorWiseSize(__FILE__, __LINE__);
//
// Utilities::PrintTuple(size1, std::cout, ", ", MpiHappyHeart().GetRankPreffix() + " Matrix1 -> [");
// Utilities::PrintTuple(size2, std::cout, ", ", MpiHappyHeart().GetRankPreffix() + " Matrix2 -> [");
// // Utilities::PrintTuple(size3, std::cout, ", ", MpiHappyHeart().GetRankPreffix() + " Matrix3 -> [");
//
// Wrappers::Petsc::MatMatMultSymbolic(interp_matrix,
// fluid_mass_matrix_on_solid_mesh,
// intermediate,
// __FILE__, __LINE__);
//
// Wrappers::Petsc::MatMatMultNumeric(interp_matrix,
// fluid_mass_matrix_on_solid_mesh,
// intermediate,
// __FILE__, __LINE__);
//
// Wrappers::Petsc::MatTransposeMatMultSymbolic(interp_matrix,
// intermediate,
// t11,
// __FILE__, __LINE__);
} }
......
...@@ -342,16 +342,23 @@ namespace HappyHeart ...@@ -342,16 +342,23 @@ namespace HappyHeart
/*! /*!
* \brief Performs Matrix-Matrix Multiplication C = A^T * B. * \class doxygen_hide_mat_transpose_mat_mult
* *
* \copydetails doxygen_hide_matmatmult_warning * \copydetails doxygen_hide_matmatmult_warning
*
* \param[in] matrix1 A in C = A^T * B. * \param[in] matrix1 A in C = A^T * B.
* \param[in] matrix2 B in C = A^T * B. * \param[in] matrix2 B in C = A^T * B.
* \param[out] matrix3 C in B in C = A^T * B. The matrix must be not allocated when this function is called. * \param[out] matrix3 C in B in C = A^T * B. The matrix must be not allocated when this function is called.
* \param[in] invoking_file File that invoked the function or class; usually __FILE__. * \param[in] invoking_file File that invoked the function or class; usually __FILE__.
* \param[in] invoking_line File that invoked the function or class; usually __LINE__. * \param[in] invoking_line File that invoked the function or class; usually __LINE__.
*/ */
/*!
* \brief Performs Matrix-Matrix Multiplication C = A^T * B.
*
* \copydetails doxygen_hide_mat_transpose_mat_mult
*
*/
template template
< <
class MatrixT, class MatrixT,
...@@ -369,6 +376,62 @@ namespace HappyHeart ...@@ -369,6 +376,62 @@ namespace HappyHeart
const char* invoking_file, int invoking_line); const char* invoking_file, int invoking_line);
/*!
* \class doxygen_hide_matrix_transpose_matrix_mult_hack
*
* \warning Petsc doesn't propose this handy function, so I hacked it with MatReuse = MAT_INITIAL_MATRIX
* in first (Symbolic) call and MatReuse = MAT_REUSE_MATRIX in subsequent calls. This works in my case...
* but is still clear ly hack.
*/
/*!
* \brief Performs Matrix-Matrix Multiplication C = A^T * B in order to compute structure of C
*
* \copydetails doxygen_hide_matrix_transpose_matrix_mult_hack
* \copydetails doxygen_hide_mat_transpose_mat_mult
*/
template
<
class MatrixT,
class MatrixU
>
std::enable_if_t
<
std::is_base_of<Private::BaseMatrix, MatrixT>::value
&& std::is_base_of<Private::BaseMatrix, MatrixU>::value,
void
>
MatTransposeMatMultSymbolic(const MatrixT& matrix1,
const MatrixT& matrix2,
MatrixU& matrix3,
const char* invoking_file, int invoking_line);
/*!
* \brief Performs Matrix-Matrix Multiplication C = A^T * B assuming C pattern is already in place.
*
* \copydetails doxygen_hide_matrix_transpose_matrix_mult_hack
* \copydetails doxygen_hide_mat_transpose_mat_mult
*/
template
<
class MatrixT,
class MatrixU
>
std::enable_if_t
<
std::is_base_of<Private::BaseMatrix, MatrixT>::value
&& std::is_base_of<Private::BaseMatrix, MatrixU>::value,
void
>
MatTransposeMatMultNumeric(const MatrixT& matrix1,
const MatrixT& matrix2,
MatrixU& matrix3,
const char* invoking_file, int invoking_line);
/*! /*!
* \brief Performs Matrix-Matrix Multiplication C = A * B^T. * \brief Performs Matrix-Matrix Multiplication C = A * B^T.
* *
......
...@@ -342,6 +342,59 @@ namespace HappyHeart ...@@ -342,6 +342,59 @@ namespace HappyHeart
} }
template
<
class MatrixT,
class MatrixU
>
inline std::enable_if_t
<
std::is_base_of<Private::BaseMatrix, MatrixT>::value
&& std::is_base_of<Private::BaseMatrix, MatrixU>::value,
void
>
MatTransposeMatMultSymbolic(const MatrixT& matrix1,
const MatrixT& matrix2,
MatrixU& matrix3,
const char* invoking_file, int invoking_line)
{
MatTransposeMatMult(matrix1,
matrix2,
matrix3,
invoking_file, invoking_line);
}
template
<
class MatrixT,
class MatrixU
>
std::enable_if_t
<
std::is_base_of<Private::BaseMatrix, MatrixT>::value
&& std::is_base_of<Private::BaseMatrix, MatrixU>::value,
void
>
MatTransposeMatMultNumeric(const MatrixT& matrix1,
const MatrixT& matrix2,
MatrixU& matrix3,
const char* invoking_file, int invoking_line)
{
Mat result = matrix3.Internal();
int error_code = ::MatTransposeMatMult(matrix1.Internal(),
matrix2.Internal(),
MAT_REUSE_MATRIX,
PETSC_DEFAULT,
&result);
if (error_code)
throw ExceptionNS::Exception(error_code, "MatTransposeMatMult", invoking_file, invoking_line);
}
/*! /*!
* \brief Performs Matrix-Matrix Multiplication C = A * B^T. * \brief Performs Matrix-Matrix Multiplication C = A * B^T.
* *
......
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