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

#1022 Now we're considering a 2x2 rather than a 3x3 matrix. However new...

#1022 Now we're considering a 2x2 rather than a 3x3 matrix. However new contribution named T21 (T31 with my not yet updated numbering) is not implemented; as a result Newton FP converges with difficulty.
parent 4d671de6
......@@ -98,7 +98,7 @@ NumberingSubset15 = {
-- Name of the numbering subset (not really used; at the moment I just need one input parameter to ground
-- the possible values to choose elsewhere).
-- Expected format: "VALUE"
name = 'fluid_mass_velocity_pressure',
name = 'fluid_monolithic',
-- Whether a vector defined on this numbering subset might be used to compute a movemesh. If true, a
-- FEltSpace featuring this numbering subset will compute additional quantities to enable fast computation.
......@@ -1066,16 +1066,16 @@ FiniteElementSpace12 = {
-- List of all unknowns defined in the finite element space. Unknowns here must be defined in this file as
-- an 'Unknown' block; expected name/identifier is the name given there.
-- Expected format: {"VALUE1", "VALUE2", ...}
unknown_list = { 'fluid_mass', 'fluid_pressure', 'fluid_velocity', 'porosity', 'solid_displacement' },
unknown_list = { 'fluid_mass', 'fluid_velocity', 'porosity', 'solid_displacement' },
-- List of the shape function to use for each unknown;
-- Expected format: {"VALUE1", "VALUE2", ...}
-- shape_function_list = { 'P2', 'P1' },
shape_function_list = { 'P1', 'P1', 'P1b', 'P1', 'P1' },
shape_function_list = { 'P1', 'P1b', 'P1', 'P1' },
-- List of the numbering subset to use for each unknown;
-- Expected format: {VALUE1, VALUE2, ...}
numbering_subset_list = { 15, 15, 15, 13, 12 }
numbering_subset_list = { 15, 15, 13, 12 }
}
-- FiniteElementSpace13 -- Fluid robin interface
......
......@@ -97,7 +97,7 @@ NumberingSubset15 = {
-- Name of the numbering subset (not really used; at the moment I just need one input parameter to ground
-- the possible values to choose elsewhere).
-- Expected format: "VALUE"
name = 'fluid_mass_velocity_pressure',
name = 'fluid_monolithic',
-- Whether a vector defined on this numbering subset might be used to compute a movemesh. If true, a
-- FEltSpace featuring this numbering subset will compute additional quantities to enable fast computation.
......
......@@ -36,7 +36,7 @@ namespace HappyHeart
decltype(auto) fluid_god_of_dof = god_of_dof_manager.GetGodOfDof(EnumUnderlyingType(MeshIndex::fluid));
decltype(auto) numbering_subset =
fluid_god_of_dof.GetNumberingSubset(EnumUnderlyingType(NumberingSubsetIndex::fluid_mass_velocity_pressure));
fluid_god_of_dof.GetNumberingSubset(EnumUnderlyingType(NumberingSubsetIndex::fluid_monolithic));
new_ = std::make_unique<GlobalVector>(numbering_subset);
auto& new_value = *new_;
......
......@@ -167,7 +167,7 @@ namespace HappyHeart
private:
///! \name Vectors defined on \a NumberingSubset NumberingSubsetIndex::fluid_mass_velocity_pressure.
///! \name Vectors defined on \a NumberingSubset NumberingSubsetIndex::fluid_monolithic.
///@{
/*!
......
......@@ -484,12 +484,8 @@ namespace HappyHeart
*/
GlobalMatrix& GetNonCstT11InMonolithic() noexcept;
GlobalMatrix& GetNonCstT21InMonolithic() noexcept;
GlobalMatrix& GetNonCstT33InMonolithic() noexcept;
GlobalMatrix& GetNonCstT22InMonolithic() noexcept;
/*!
* \brief Work matrix that gets the exact same structure as system matrix.
*
......@@ -510,14 +506,8 @@ namespace HappyHeart
* pattern (an assert is nonetheless there to check it does respect the pattern delimited by the
* numbering subsets).
*/
GlobalMatrix& GetNonCstT32InMonolithic() noexcept;
GlobalMatrix& GetNonCstIntermediateT32() noexcept;
GlobalMatrix& GetNonCstT32Block() noexcept;
GlobalMatrix& GetNonCstT33Block() noexcept;
GlobalMatrix& GetNonCstT22Block() noexcept;
/*!
* \brief Non constant accessor to the matrix into which mass operator for fluid mass / fluid pressure
......@@ -526,16 +516,10 @@ namespace HappyHeart
*/
GlobalMatrix& GetNonCstFluidMassPressureMatrixOnSolidMesh() noexcept;
//! Non constant accessor to the t21 on the non monolithic FEltSpace.
GlobalMatrix& GetNonCstT21OnFluidMesh() noexcept;
///@}
const ::HappyHeart::GlobalVariationalOperatorNS::TransientSource<ParameterNS::Type::scalar>&
GetFluidPressureSolidMeshOperator() const noexcept;
public:
/*!
......@@ -576,26 +560,15 @@ namespace HappyHeart
void ComputeT11(DoRecomputeBlock do_recompute_block,
Wrappers::Petsc::DoReuseMatrix do_reuse_matrix = Wrappers::Petsc::DoReuseMatrix::yes);
//! Compute T21.
void ComputeT21(DoRecomputeBlock do_recompute_block,
Wrappers::Petsc::DoReuseMatrix do_reuse_matrix = Wrappers::Petsc::DoReuseMatrix::yes);
//! Compute T13.
void ComputeT13(DoRecomputeBlock do_recompute_block,
Wrappers::Petsc::DoReuseMatrix do_reuse_matrix = Wrappers::Petsc::DoReuseMatrix::yes);
//! Compute T32.
void ComputeT32(DoRecomputeBlock do_recompute_block,
Wrappers::Petsc::DoReuseMatrix do_reuse_matrix = Wrappers::Petsc::DoReuseMatrix::yes);
//! Compute T33.
void ComputeT33(DoApplyBoundaryCondition do_apply_boundary_condition,
DoRecomputeBlock do_recompute_block);
//! Compute T22.
void ComputeT22(DoRecomputeBlock do_recompute_block);
private:
......@@ -681,15 +654,6 @@ namespace HappyHeart
*/
GlobalMatrix::unique_ptr t11_in_monolithic_ = nullptr;
/*!
* \brief Work matrix that gets the exact same structure as system matrix.
*
* This matrix isn't allocated directly: the MatMatMult functions are actually writing its
* pattern (an assert is nonetheless there to check it does respect the pattern delimited by the
* numbering subsets).
*/
GlobalMatrix::unique_ptr t21_in_monolithic_ = nullptr;
/*!
* \brief Work matrix that gets the exact same structure as system matrix.
*
......@@ -704,40 +668,12 @@ namespace HappyHeart
GlobalMatrix::unique_ptr t13_block_ = nullptr;
/*!
* \brief Work matrix that gets the exact same structure as system matrix.
*
* This matrix isn't allocated directly: the MatMatMult functions are actually writing its
* pattern (an assert is nonetheless there to check it does respect the pattern delimited by the
* numbering subsets).
*/
GlobalMatrix::unique_ptr t32_in_monolithic_ = nullptr;
GlobalMatrix::unique_ptr intermediate_t32_ = nullptr;
GlobalMatrix::unique_ptr t32_block_ = nullptr;
//! T21 on the non monolithic \a FEltSpace.
GlobalMatrix::unique_ptr t21_on_fluid_mesh_ = nullptr;
//! Intermediate work matrix to transfer T21 from solid to fluid.
GlobalMatrix::unique_ptr intermediate_t21_ = nullptr;
//! Intermediate work matrix to transfer T21 from block to monolithic.
GlobalMatrix::unique_ptr intermediate_t21_monolithic_ = nullptr;
GlobalMatrix::unique_ptr t33_in_monolithic_ = nullptr;
// \todo #820 TMP for debug.
GlobalMatrix::unique_ptr t33_block_ = nullptr;
GlobalMatrix::unique_ptr t23_block_ = nullptr;
GlobalMatrix::unique_ptr t22_block_ = nullptr;
GlobalMatrix::unique_ptr t22_in_monolithic_ = nullptr;
///@}
/*!
......
......@@ -63,7 +63,7 @@ namespace HappyHeart
decltype(auto) god_of_dof = parent::GetGodOfDof();
decltype(auto) numbering_subset =
god_of_dof.GetNumberingSubset(EnumUnderlyingType(NumberingSubsetIndex::fluid_mass_velocity_pressure));
god_of_dof.GetNumberingSubset(EnumUnderlyingType(NumberingSubsetIndex::fluid_monolithic));
parent::SolveNonLinear(numbering_subset, numbering_subset,
__FILE__, __LINE__);
......@@ -134,7 +134,7 @@ namespace HappyHeart
monolithic_data.SetNew(corrected_values);
decltype(auto) god_of_dof = parent::GetGodOfDof();
decltype(auto) numbering_subset =
god_of_dof.GetNumberingSubset(EnumUnderlyingType(NumberingSubsetIndex::fluid_mass_velocity_pressure));
god_of_dof.GetNumberingSubset(EnumUnderlyingType(NumberingSubsetIndex::fluid_monolithic));
parent::GetNonCstSystemSolution(numbering_subset).Copy(corrected_values, __FILE__, __LINE__);
......@@ -168,7 +168,7 @@ namespace HappyHeart
{
decltype(auto) god_of_dof = parent::GetGodOfDof();
decltype(auto) numbering_subset =
god_of_dof.GetNumberingSubset(EnumUnderlyingType(NumberingSubsetIndex::fluid_mass_velocity_pressure));
god_of_dof.GetNumberingSubset(EnumUnderlyingType(NumberingSubsetIndex::fluid_monolithic));
decltype(auto) variable_holder = this->GetNonCstVariableHolder();
......@@ -210,7 +210,7 @@ namespace HappyHeart
decltype(auto) god_of_dof = parent::GetGodOfDof();
decltype(auto) numbering_subset =
god_of_dof.GetNumberingSubset(EnumUnderlyingType(NumberingSubsetIndex::fluid_mass_velocity_pressure));
god_of_dof.GetNumberingSubset(EnumUnderlyingType(NumberingSubsetIndex::fluid_monolithic));
decltype(auto) system_matrix = parent::GetSystemMatrix(numbering_subset, numbering_subset);
......@@ -249,7 +249,7 @@ namespace HappyHeart
{
decltype(auto) god_of_dof = parent::GetGodOfDof();
decltype(auto) numbering_subset =
god_of_dof.GetNumberingSubset(EnumUnderlyingType(NumberingSubsetIndex::fluid_mass_velocity_pressure));
god_of_dof.GetNumberingSubset(EnumUnderlyingType(NumberingSubsetIndex::fluid_monolithic));
auto& system_matrix = parent::GetNonCstSystemMatrix(numbering_subset, numbering_subset);
system_matrix.ZeroEntries(__FILE__, __LINE__);
......@@ -262,11 +262,8 @@ namespace HappyHeart
// IMPORTANT: T33 must be assembled first, as there is a Dirichlet boundary condition that deals only
// with its dofs.
ComputeT33(DoApplyBoundaryCondition::yes, DoRecomputeBlock::yes);
ComputeT22(DoRecomputeBlock::yes);
ComputeT13(DoRecomputeBlock::yes);
ComputeT11(DoRecomputeBlock::yes);
ComputeT21(DoRecomputeBlock::yes);
ComputeT32(DoRecomputeBlock::yes);
break;
}
......@@ -275,11 +272,8 @@ namespace HappyHeart
// IMPORTANT: T33 must be assembled first, as there is a Dirichlet boundary condition that deals only
// with its dofs.
ComputeT33(DoApplyBoundaryCondition::yes, DoRecomputeBlock::yes);
ComputeT22(DoRecomputeBlock::yes);
ComputeT13(DoRecomputeBlock::yes);
ComputeT11(DoRecomputeBlock::no);
ComputeT21(DoRecomputeBlock::yes);
ComputeT32(DoRecomputeBlock::yes);
break;
}
......@@ -288,22 +282,16 @@ namespace HappyHeart
// IMPORTANT: T33 must be assembled first, as there is a Dirichlet boundary condition that deals only
// with its dofs.
ComputeT33(DoApplyBoundaryCondition::yes, DoRecomputeBlock::no);
ComputeT22(DoRecomputeBlock::no);
ComputeT13(DoRecomputeBlock::no);
ComputeT11(DoRecomputeBlock::no);
ComputeT21(DoRecomputeBlock::yes);
ComputeT32(DoRecomputeBlock::no);
break;
}
case BlocksToRecompute::t33_without_dirichlet:
{
ComputeT33(DoApplyBoundaryCondition::no, DoRecomputeBlock::yes);
ComputeT22(DoRecomputeBlock::no);
ComputeT13(DoRecomputeBlock::no);
ComputeT11(DoRecomputeBlock::no);
ComputeT21(DoRecomputeBlock::no);
ComputeT32(DoRecomputeBlock::no);
break;
}
......@@ -347,7 +335,7 @@ namespace HappyHeart
decltype(auto) god_of_dof = parent::GetGodOfDof();
decltype(auto) bc_manager = DirichletBoundaryConditionManager::GetInstance();
decltype(auto) numbering_subset =
god_of_dof.GetNumberingSubset(EnumUnderlyingType(NumberingSubsetIndex::fluid_mass_velocity_pressure));
god_of_dof.GetNumberingSubset(EnumUnderlyingType(NumberingSubsetIndex::fluid_monolithic));
auto& rhs = parent::GetNonCstSystemRhs(numbering_subset);
rhs.ZeroEntries(__FILE__, __LINE__);
......@@ -370,73 +358,7 @@ namespace HappyHeart
variable_holder.template DevPrint<DevPhase::newton_fp>(rhs,
"darcy");
// \todo #820 Add the contribution from pressure to the rhs.
{
// First the dWdm contribution:
typename DataNS::NewFluidPressure<HyperelasticLawT>::global_vector_temporary_solid_manager
pressure_source_on_solid_helper(fluid_pressure_data);
auto& pressure_source_on_solid = pressure_source_on_solid_helper.GetNonCstVector();
pressure_source_on_solid.ZeroEntries(__FILE__, __LINE__);
GlobalVectorWithCoefficient vector_with_coeff(pressure_source_on_solid, 1.);
// Then add the current pressure on solid contribution.
typename DataNS::NewFluidPressure<HyperelasticLawT>::global_vector_temporary_solid_manager
pressure_source_mass_term_helper(fluid_pressure_data);
auto& pressure_source_mass_term = pressure_source_mass_term_helper.GetNonCstVector();
pressure_source_mass_term.ZeroEntries(__FILE__, __LINE__);
decltype(auto) time_manager = parent::GetTimeManager();
decltype(auto) time = time_manager.GetTime();
GlobalVectorWithCoefficient fluid_pressure_solid_mesh_with_coeff(pressure_source_mass_term, 1.);
fluid_pressure_solid_mesh_operator_->Assemble(std::make_tuple(std::ref(fluid_pressure_solid_mesh_with_coeff)),
time);
variable_holder.template DevPrint<DevPhase::newton_fp>(pressure_source_mass_term,
"pressure_contrib_to_rhs_only_source_2");
Wrappers::Petsc::AXPY(1.,
pressure_source_mass_term,
pressure_source_on_solid,
__FILE__, __LINE__);
variable_holder.template DevPrint<DevPhase::newton_fp>(pressure_source_on_solid,
"pressure_contrib_to_rhs");
// Project onto monolithic and add contribution to rhs.
DataNS::Monolithic::global_vector_temporary_manager work(monolithic_data);
auto& monolithic_pressure_source = work.GetNonCstVector();
Wrappers::Petsc::MatMult(interpolator_holder.GetMatrixFluidPressureOnSolidToMonolithic(),
pressure_source_on_solid,
monolithic_pressure_source,
__FILE__, __LINE__);
variable_holder.template DevPrint<DevPhase::newton_fp>(pressure_source_on_solid,
"monolithic_pressure_contrib_to_rhs");
Wrappers::Petsc::AXPY(1.,
monolithic_pressure_source,
rhs,
__FILE__, __LINE__);
// variable_holder.template DevPrint<DevPhase::newton_fp>(fluid_pressure_data.GetOnSolid(),
// "pressure_on_solid");
variable_holder.template DevPrint<DevPhase::newton_fp>(solid_displacement_data.GetCurrent(),
"newtonFP_dispSr_old");
variable_holder.template DevPrint<DevPhase::newton_fp>(fluid_mass_data.GetNew(),
"newtonFP_fluidmass");
variable_holder.template DevPrint<DevPhase::newton_fp>(fluid_mass_data.GetCurrent(),
"newtonFP_fluidmass_old");
}
// \todo #1022 Check in Bruno's note if contribution is to be added here.
{
DataNS::Monolithic::global_vector_temporary_manager work(monolithic_data);
auto& contribution = work.GetNonCstVector();
......
......@@ -84,17 +84,7 @@ namespace HappyHeart
return *mass_operator_fluid_mass_on_solid_mesh_;
}
template<class HyperelasticLawT>
inline const ::HappyHeart::GlobalVariationalOperatorNS::TransientSource<ParameterNS::Type::scalar>&
VariationalFormulation<HyperelasticLawT>
::GetFluidPressureOnSolidMeshOperator() const noexcept
{
assert(!(!fluid_pressure_solid_mesh_operator_));
return *fluid_pressure_solid_mesh_operator_;
}
template<class HyperelasticLawT>
inline const GlobalMatrix& VariationalFormulation<HyperelasticLawT>
::GetFluidMassMatrixOnSolidMesh() const noexcept
......@@ -148,15 +138,6 @@ namespace HappyHeart
}
template<class HyperelasticLawT>
inline GlobalMatrix& VariationalFormulation<HyperelasticLawT>
::GetNonCstT21InMonolithic() noexcept
{
assert(!(!t21_in_monolithic_));
return *t21_in_monolithic_;
}
template<class HyperelasticLawT>
inline GlobalMatrix& VariationalFormulation<HyperelasticLawT>
::GetNonCstT13InMonolithic() noexcept
......@@ -184,41 +165,6 @@ namespace HappyHeart
}
template<class HyperelasticLawT>
inline GlobalMatrix& VariationalFormulation<HyperelasticLawT>
::GetNonCstT32InMonolithic() noexcept
{
assert(!(!t32_in_monolithic_));
return *t32_in_monolithic_;
}
template<class HyperelasticLawT>
inline GlobalMatrix& VariationalFormulation<HyperelasticLawT>
::GetNonCstT22InMonolithic() noexcept
{
assert(!(!t22_in_monolithic_));
return *t22_in_monolithic_;
}
template<class HyperelasticLawT>
inline GlobalMatrix& VariationalFormulation<HyperelasticLawT>
::GetNonCstIntermediateT32() noexcept
{
assert(!(!intermediate_t32_));
return *intermediate_t32_;
}
template<class HyperelasticLawT>
inline GlobalMatrix& VariationalFormulation<HyperelasticLawT>::GetNonCstT32Block() noexcept
{
assert(!(!t32_block_));
return *t32_block_;
}
template<class HyperelasticLawT>
inline GlobalMatrix& VariationalFormulation<HyperelasticLawT>::GetNonCstT33Block() noexcept
......@@ -228,13 +174,6 @@ namespace HappyHeart
}
template<class HyperelasticLawT>
inline GlobalMatrix& VariationalFormulation<HyperelasticLawT>::GetNonCstT22Block() noexcept
{
assert(!(!t22_block_));
return *t22_block_;
}
template<class HyperelasticLawT>
inline const ::HappyHeart::GlobalVariationalOperatorNS::Mass& VariationalFormulation<HyperelasticLawT>
......@@ -270,33 +209,14 @@ namespace HappyHeart
return *fluid_mass_pressure_matrix_on_solid_mesh_;
}
template<class HyperelasticLawT>
inline const GlobalVariationalOperatorNS::T21<HyperelasticLawT>& VariationalFormulation<HyperelasticLawT>
::GetT21Operator() const noexcept
{
assert(!(!t21_operator_));
return *t21_operator_;
}
template<class HyperelasticLawT>
inline const HyperelasticLawT& VariationalFormulation<HyperelasticLawT>
::GetHyperelasticLaw() const noexcept
{
return hyperelastic_law_;
}
template<class HyperelasticLawT>
inline GlobalMatrix& VariationalFormulation<HyperelasticLawT>::GetNonCstT21OnFluidMesh() noexcept
{
assert(!(!t21_on_fluid_mesh_));
return *t21_on_fluid_mesh_;
}
template<class HyperelasticLawT>
inline const GlobalVariationalOperatorNS::HybridVector<HyperelasticLawT>&
......@@ -374,7 +294,7 @@ namespace HappyHeart
template<class HyperelasticLawT>
inline const ::HappyHeart::GlobalVariationalOperatorNS::TransientSource<ParameterNS::Type::scalar>&
VariationalFormulation<HyperelasticLawT>
::GetFluidPressureSolidMeshOperator() const noexcept
::GetFluidPressureOnSolidMeshOperator() const noexcept
{
decltype(auto) new_fluid_pressure_data = this->GetNonCstNewFluidPressureData();
......
......@@ -39,7 +39,7 @@ namespace HappyHeart
decltype(auto) fluid_god_of_dof = parent::GetGodOfDof();
decltype(auto) numbering_subset =
fluid_god_of_dof.GetNumberingSubset(EnumUnderlyingType(NumberingSubsetIndex::fluid_mass_velocity_pressure));
fluid_god_of_dof.GetNumberingSubset(EnumUnderlyingType(NumberingSubsetIndex::fluid_monolithic));
auto& system_matrix = parent::GetNonCstSystemMatrix(numbering_subset, numbering_subset);
decltype(auto) interpolator_holder = this->GetInterpolatorHolder();
......@@ -128,99 +128,6 @@ namespace HappyHeart
}
template<class HyperelasticLawT>
void VariationalFormulation<HyperelasticLawT>
::ComputeT21(DoRecomputeBlock do_recompute_block,
Wrappers::Petsc::DoReuseMatrix do_reuse_matrix)
{
// \todo #820 Currently two transferts are done for dev purposes; however it is possible to reduce
// that to only one!
decltype(auto) god_of_dof = parent::GetGodOfDof();
decltype(auto) numbering_subset =
god_of_dof.GetNumberingSubset(EnumUnderlyingType(NumberingSubsetIndex::fluid_mass_velocity_pressure));
auto& t21_in_monolithic = *t21_in_monolithic_;
auto& system_matrix = parent::GetNonCstSystemMatrix(numbering_subset, numbering_subset);
decltype(auto) interpolator_holder = this->GetInterpolatorHolder();
switch(do_recompute_block)
{
case DoRecomputeBlock::no:
break;
case DoRecomputeBlock::yes:
{
auto& t21_on_solid_mesh = GetNonCstFluidMassPressureMatrixOnSolidMesh();
t21_on_solid_mesh.ZeroEntries(__FILE__, __LINE__);
{
GlobalMatrixWithCoefficient matrix_with_coeff(t21_on_solid_mesh, 1.);
GetT21Operator().Assemble(std::make_tuple(std::ref(matrix_with_coeff)));
}
auto& t21_on_fluid_mesh = GetNonCstT21OnFluidMesh();