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

#723 FSI/Newton: tyring to use Petsc MatShell feature. Neither Mumps nor...

#723 FSI/Newton: tyring to use Petsc MatShell feature. Neither Mumps nor UMFPack are actually supported for it; I'll try to use GMRES.
parent 29075d9a
......@@ -626,6 +626,45 @@ Petsc1 = {
} -- Petsc
Petsc2 = {
-- Absolute tolerance
-- Expected format: {VALUE1, VALUE2, ...}
-- Constraint: v > 0.
absoluteTolerance = 1e-10,
-- gmresStart
-- Expected format: {VALUE1, VALUE2, ...}
-- Constraint: v >= 0
gmresRestart = 200,
-- Maximum iteration
-- Expected format: {VALUE1, VALUE2, ...}
-- Constraint: v > 0
maxIteration = 1000,
-- List of preconditioner: { none jacobi sor lu bjacobi ilu asm cholesky }.
-- To use mumps:
-- preconditioner = lu
-- Expected format: {"VALUE1", "VALUE2", ...}
-- Constraint: ops_in(v, 'lu')
preconditioner = 'lu',
-- Relative tolerance
-- Expected format: {VALUE1, VALUE2, ...}
-- Constraint: v > 0.
relativeTolerance = 1e-6,
-- List of solver: { chebychev cg gmres preonly bicg python };
-- Expected format: {"VALUE1", "VALUE2", ...}
-- Constraint: ops_in(v, {'Mumps', 'Umfpack'})
solver = 'Mumps',
} -- Petsc
Fluid = {
Viscosity = {
......
......@@ -87,8 +87,8 @@ namespace HappyHeart
template<IsFactorized IsFactorizedT>
void VariationalFormulation<DerivedT, SolverIndexT>
::SolveLinear(const GlobalMatrix& matrix,
const GlobalVector& rhs,
GlobalVector& out)
const GlobalVector& rhs,
GlobalVector& out)
{
assert(out.GetNumberingSubset() == rhs.GetNumberingSubset());
......
......@@ -71,7 +71,8 @@ namespace HappyHeart
enum class SolverIndex
{
main_solver = 1
main_solver = 1,
dh_solver
};
......@@ -140,7 +141,8 @@ namespace HappyHeart
InputParameter::FEltSpace<EnumUnderlyingType(FEltSpaceIndex::inlet_border)>,
InputParameter::Petsc<EnumUnderlyingType(SolverIndex::main_solver)>,
InputParameter::Petsc<EnumUnderlyingType(SolverIndex::dh_solver)>,
InputParameter::Fluid::Viscosity,
InputParameter::Fluid::Density,
......
......@@ -185,6 +185,8 @@ namespace HappyHeart
///@}
public: // \todo #723 See what is better once shell has been mastered.
void dH();
......@@ -251,11 +253,11 @@ namespace HappyHeart
//! Solid velocity on interface.
GlobalVector& GetNonCstSolidVelocityOnInterface() noexcept;
//! Access to the linear solver used to solve the diffferential of the solid variational formulation.
Wrappers::Petsc::Snes& GetNonCstDifferentialSolidSolver() noexcept;
Wrappers::Petsc::Snes& GetNonCstdHSolver() noexcept;
private:
......@@ -305,6 +307,8 @@ namespace HappyHeart
ConformInterpolatorNS::P1_to_P2::unique_ptr solid_to_fluid_velocity_on_fsi_interpolator_ = nullptr;
Wrappers::Petsc::Snes::unique_ptr differential_solid_solver_ = nullptr;
Wrappers::Petsc::Snes::unique_ptr dH_solver_ = nullptr;
};
......
......@@ -246,6 +246,7 @@ namespace HappyHeart
}
differential_solid_solver_ = InitSolver<EnumUnderlyingType(SolverIndex::main_solver)>(mpi, input_parameter_data);
dH_solver_ = InitSolver<EnumUnderlyingType(SolverIndex::dh_solver)>(mpi, input_parameter_data);
}
......
......@@ -168,6 +168,14 @@ namespace HappyHeart
}
template<class SolidVariationalFormulationPolicyT>
Wrappers::Petsc::Snes& Model<SolidVariationalFormulationPolicyT>::GetNonCstdHSolver() noexcept
{
assert(!(!dH_solver_));
return *dH_solver_;
}
template<class SolidVariationalFormulationPolicyT>
const GlobalVector& Model<SolidVariationalFormulationPolicyT>::GetSolidVelocityOnInterface() const noexcept
{
......
......@@ -18,6 +18,21 @@ namespace HappyHeart
namespace FSI_EINS
{
template<class SolidVariationalFormulationPolicyT>
PetscInt ApplydH(Mat shell_matrix, Vec, Vec)
{
Model<SolidVariationalFormulationPolicyT>* context;
int error_code = MatShellGetContext(shell_matrix, &context);
context->dH();
return 0;
}
template<class SolidVariationalFormulationPolicyT>
......@@ -48,13 +63,24 @@ namespace HappyHeart
UpdateSolidDisplacementAfterSoF();
Wrappers::Petsc::ShellMatrix<std::false_type, self>(mpi,
0u, 0u,
0u, 0u,
this,
__FILE__, __LINE__);
dH();
const auto processor_wise_size = GetSolidDisplacementAfterSoF().GetProcessorWiseSize(__FILE__, __LINE__);
const auto program_wise_size = GetSolidDisplacementAfterSoF().GetProgramWiseSize(__FILE__, __LINE__);
Wrappers::Petsc::ShellMatrix<std::false_type, self> test(mpi,
processor_wise_size, processor_wise_size,
program_wise_size, program_wise_size,
this,
__FILE__, __LINE__);
auto petsc_mat = test.Internal();
int error_code = MatShellSetOperation(petsc_mat, MATOP_MULT,
(void(*)(void))ApplydH<SolidVariationalFormulationPolicyT>);
dH_solver_->SolveLinear(test, test,
GetNonCstSolidDisplacementAfterSoF(),
GetNonCstSolidDisplacementAfterSoF(),
__FILE__, __LINE__);
}
......
......@@ -205,120 +205,7 @@ namespace HappyHeart
return preconditioner;
}
void Snes::SolveLinear(const Matrix& matrix,
const Matrix& preconditioner,
const Vector& rhs,
Vector& solution,
const char* invoking_file, int invoking_line)
{
KSP ksp = GetKsp(invoking_file, invoking_line);
TimeKeep::GetInstance().PrintTimeElapsed("BeforeSetOperator");
int error_code = KSPSetOperators(ksp, matrix.Internal(), preconditioner.Internal());
if (error_code)
throw ExceptionNS::Exception(error_code, "KSPSetOperators", invoking_file, invoking_line);
// If there was a previous call to this very function, this argument is set to PETSC_TRUE.
// That's not what we want: matrix and preconditioner may not be the same as in previous call.
error_code = KSPSetReusePreconditioner(ksp, PETSC_FALSE);
if (error_code)
throw ExceptionNS::Exception(error_code, "KSPSetReusePreconditioner", invoking_file, invoking_line);
TimeKeep::GetInstance().PrintTimeElapsed("AfterSetOperator");
auto pc = GetPreconditioner(__FILE__, __LINE__);
error_code = PCFactorSetUpMatSolverPackage(pc); /* call MatGetFactor() to create F */
if (error_code)
throw ExceptionNS::Exception(error_code, "PCFactorSetUpMatSolverPackage",
invoking_file, invoking_line);
Mat F;
error_code = PCFactorGetMatrix(pc, &F);
if (error_code)
throw ExceptionNS::Exception(error_code, "PCFactorGetMatrix",
invoking_file, invoking_line);
// /* sequential ordering */
PetscInt icntl = 7, ival = 2;
PetscScalar val;
error_code = MatMumpsSetIcntl(F, icntl, ival);
if (error_code)
throw ExceptionNS::Exception(error_code, "MatMumpsSetIcntl",
invoking_file, invoking_line);
/* threshhold for row pivot detection */
error_code = MatMumpsSetIcntl(F,24,1);
if (error_code)
throw ExceptionNS::Exception(error_code, "MatMumpsSetIcntl",
invoking_file, invoking_line);
icntl = 3; val = 1.e-6;
error_code = MatMumpsSetCntl(F,icntl,val);
if (error_code)
throw ExceptionNS::Exception(error_code, "MatMumpsSetIcntl",
invoking_file, invoking_line);
/* compute determinant of A */
error_code = MatMumpsSetIcntl(F,33,1);
if (error_code)
throw ExceptionNS::Exception(error_code, "MatMumpsSetIcntl",
invoking_file, invoking_line);
TimeKeep::GetInstance().PrintTimeElapsed("AfterMUMPSSet");
// This optional line (would be called in KSPSolve) can lead to more explicit messages.
// Petsc documentation quote: 'The explicit call of this routine enables the separate monitoring of any
// computations performed during the set up phase, such as incomplete factorization for the ILU
// preconditioner.'
error_code = ::KSPSetUp(ksp);
if (error_code)
throw ExceptionNS::Exception(error_code, "KSPSetUp", invoking_file, invoking_line);
// It seems I need this step to ensure Petsc doesn't do again the preconditioning; this is weird
// because Petsc documentation doesn't mention it at all (this function is there in cases we want
// to force to reuse the preconditioning even if the matrix or preconditioner has changed).
error_code = KSPSetReusePreconditioner(ksp, PETSC_TRUE);
if (error_code)
throw ExceptionNS::Exception(error_code, "KSPSetReusePreconditioner", invoking_file, invoking_line);
TimeKeep::GetInstance().PrintTimeElapsed("AfterSetUp");
error_code = ::KSPSolve(ksp, rhs.Internal(), solution.Internal());
if (error_code)
throw ExceptionNS::Exception(error_code, "KSPSolve", invoking_file, invoking_line);
TimeKeep::GetInstance().PrintTimeElapsed("AfterKSPSolve");
PetscInt Niteration;
error_code = KSPGetIterationNumber(ksp, &Niteration);
if (error_code)
throw ExceptionNS::Exception(error_code, "KSPGetIterationNumber", invoking_file, invoking_line);
PetscReal residual;
error_code = KSPGetResidualNorm(ksp, &residual);
if (error_code)
throw ExceptionNS::Exception(error_code, "KSPGetResidualNorm", invoking_file, invoking_line);
PrintMessageOnFirstProcessor("%s-%s converged in %d iterations (residual = %e)\n",
GetMpi(),
invoking_file, invoking_line,
GetKspType(invoking_file, invoking_line).c_str(),
GetPreconditionerType(invoking_file, invoking_line).c_str(),
static_cast<int>(Niteration),
static_cast<double>(residual));
}
void Snes::SolveLinear(const Vector& rhs,
Vector& solution,
......
......@@ -17,6 +17,12 @@
# include "ThirdParty/Wrappers/Petsc/Vector/Vector.hpp"
# include "ThirdParty/Wrappers/Mpi/Mpi.hpp"
#include "ThirdParty/Wrappers/Petsc/Snes.hpp"
#include "ThirdParty/Wrappers/Petsc/Print.hpp"
#include "ThirdParty/Wrappers/Petsc/Exceptions/Petsc.hpp"
#include "Utilities/TimeKeep/TimeKeep.hpp"
namespace HappyHeart
{
......@@ -135,9 +141,11 @@ namespace HappyHeart
* \param[in] invoking_line File that invoked the function or class; usually __LINE__.
*
*/
void SolveLinear(const Matrix& matrix, const Matrix& preconditioner,
const Vector& rhs, Vector& solution,
const char* invoking_file, int invoking_line);
template<class MatrixT>
std::enable_if_t<std::is_base_of<Private::BaseMatrix, MatrixT>::value, void>
SolveLinear(const MatrixT& matrix, const MatrixT& preconditioner,
const Vector& rhs, Vector& solution,
const char* invoking_file, int invoking_line);
/*!
......
......@@ -45,8 +45,121 @@ namespace HappyHeart
}
template<class MatrixT>
std::enable_if_t<std::is_base_of<Private::BaseMatrix, MatrixT>::value, void>
Snes::SolveLinear(const MatrixT& matrix,
const MatrixT& preconditioner,
const Vector& rhs,
Vector& solution,
const char* invoking_file, int invoking_line)
{
KSP ksp = GetKsp(invoking_file, invoking_line);
TimeKeep::GetInstance().PrintTimeElapsed("BeforeSetOperator");
int error_code = KSPSetOperators(ksp, matrix.Internal(), preconditioner.Internal());
if (error_code)
throw ExceptionNS::Exception(error_code, "KSPSetOperators", invoking_file, invoking_line);
// If there was a previous call to this very function, this argument is set to PETSC_TRUE.
// That's not what we want: matrix and preconditioner may not be the same as in previous call.
error_code = KSPSetReusePreconditioner(ksp, PETSC_FALSE);
if (error_code)
throw ExceptionNS::Exception(error_code, "KSPSetReusePreconditioner", invoking_file, invoking_line);
TimeKeep::GetInstance().PrintTimeElapsed("AfterSetOperator");
auto pc = GetPreconditioner(__FILE__, __LINE__);
error_code = PCFactorSetUpMatSolverPackage(pc); /* call MatGetFactor() to create F */
if (error_code)
throw ExceptionNS::Exception(error_code, "PCFactorSetUpMatSolverPackage",
invoking_file, invoking_line);
Mat F;
error_code = PCFactorGetMatrix(pc, &F);
if (error_code)
throw ExceptionNS::Exception(error_code, "PCFactorGetMatrix",
invoking_file, invoking_line);
// /* sequential ordering */
PetscInt icntl = 7, ival = 2;
PetscScalar val;
error_code = MatMumpsSetIcntl(F, icntl, ival);
if (error_code)
throw ExceptionNS::Exception(error_code, "MatMumpsSetIcntl",
invoking_file, invoking_line);
/* threshhold for row pivot detection */
error_code = MatMumpsSetIcntl(F,24,1);
if (error_code)
throw ExceptionNS::Exception(error_code, "MatMumpsSetIcntl",
invoking_file, invoking_line);
icntl = 3; val = 1.e-6;
error_code = MatMumpsSetCntl(F,icntl,val);
if (error_code)
throw ExceptionNS::Exception(error_code, "MatMumpsSetIcntl",
invoking_file, invoking_line);
/* compute determinant of A */
error_code = MatMumpsSetIcntl(F,33,1);
if (error_code)
throw ExceptionNS::Exception(error_code, "MatMumpsSetIcntl",
invoking_file, invoking_line);
TimeKeep::GetInstance().PrintTimeElapsed("AfterMUMPSSet");
// This optional line (would be called in KSPSolve) can lead to more explicit messages.
// Petsc documentation quote: 'The explicit call of this routine enables the separate monitoring of any
// computations performed during the set up phase, such as incomplete factorization for the ILU
// preconditioner.'
error_code = ::KSPSetUp(ksp);
if (error_code)
throw ExceptionNS::Exception(error_code, "KSPSetUp", invoking_file, invoking_line);
// It seems I need this step to ensure Petsc doesn't do again the preconditioning; this is weird
// because Petsc documentation doesn't mention it at all (this function is there in cases we want
// to force to reuse the preconditioning even if the matrix or preconditioner has changed).
error_code = KSPSetReusePreconditioner(ksp, PETSC_TRUE);
if (error_code)
throw ExceptionNS::Exception(error_code, "KSPSetReusePreconditioner", invoking_file, invoking_line);
TimeKeep::GetInstance().PrintTimeElapsed("AfterSetUp");
error_code = ::KSPSolve(ksp, rhs.Internal(), solution.Internal());
if (error_code)
throw ExceptionNS::Exception(error_code, "KSPSolve", invoking_file, invoking_line);
TimeKeep::GetInstance().PrintTimeElapsed("AfterKSPSolve");
PetscInt Niteration;
error_code = KSPGetIterationNumber(ksp, &Niteration);
if (error_code)
throw ExceptionNS::Exception(error_code, "KSPGetIterationNumber", invoking_file, invoking_line);
PetscReal residual;
error_code = KSPGetResidualNorm(ksp, &residual);
if (error_code)
throw ExceptionNS::Exception(error_code, "KSPGetResidualNorm", invoking_file, invoking_line);
PrintMessageOnFirstProcessor("%s-%s converged in %d iterations (residual = %e)\n",
GetMpi(),
invoking_file, invoking_line,
GetKspType(invoking_file, invoking_line).c_str(),
GetPreconditionerType(invoking_file, invoking_line).c_str(),
static_cast<int>(Niteration),
static_cast<double>(residual));
}
} // namespace Petsc
......
Supports Markdown
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