Skip to content

GitLab

Commit 997dc1b1 authored by GILLES Sebastien's avatar GILLES Sebastien
Browse files

#820 In progress: computing the whole Dirichlet vector before NewtonFP loop. A... 

#820 In progress: computing the whole Dirichlet vector before NewtonFP loop. A likely bug was found in velocity mesh -> robin operator (the numbering subset given was probably wrong).
parent 02c93349
Hide whitespace changes
Inline Side-by-side
Sources/ModelInstances/UnderDevelopment/Poromechanics/ImplicitStepFluid/InnerLoop/VariationalFormulation.hpp
View file @ 997dc1b1
......@@ -323,12 +323,7 @@ namespace HappyHeart
//! Non constant accessor to the \#820 [DEV] Vector defined on fluid mass numbering subset.
GlobalVector& GetNonCstWorkFluidMass() noexcept;
//! Constant accessor to the \#820 [DEV] Vector defined on fluid velocity numbering subset.
const GlobalVector& GetWorkFluidVelocity() const noexcept;
//! Non constant accessor to the \#820 [DEV] Vector defined on fluid velocity numbering subset.
GlobalVector& GetNonCstWorkFluidVelocity() noexcept;
//! Constant accessor to the mass div Velocity operator.
const GlobalVariationalOperatorNS::ScalarDivVectorial& GetMassDivVelocity() const noexcept;
......@@ -732,9 +727,6 @@ namespace HappyHeart
//! \#820 [DEV] Vector defined on fluid mass numbering subset.
GlobalVector::unique_ptr work_fluid_mass_ = nullptr;
//! \#820 [DEV] Vector defined on fluid velocity numbering subset.
GlobalVector::unique_ptr work_fluid_velocity_ = nullptr;
//! \#820 [DEV] Vector defined on the numbering subset for solid on fluid.
GlobalVector::unique_ptr work_solid_velocity_ = nullptr;
......
Sources/ModelInstances/UnderDevelopment/Poromechanics/ImplicitStepFluid/InnerLoop/VariationalFormulation.hxx
View file @ 997dc1b1
......@@ -66,6 +66,10 @@ namespace HappyHeart
corrected_values,
__FILE__, __LINE__);
auto& vel_contrib_on_vel_ns = variable_holder.GetVelSHalfVhfZeroedOutOfRobinInterface();
variable_holder.template DevPrint<DevPhase::newton_fp>(variable_holder.GetVelSHalfVhf(),
"dirichlet_just_before_newton_fp_loop");
......@@ -589,9 +593,10 @@ namespace HappyHeart
GlobalVectorWithCoefficient contribution_with_coeff(contribution, 1.);
decltype(auto) interpolator_holder = this->GetInterpolatorHolder();
decltype(auto) variable_holder = this->GetVariableHolder();
{
auto& vel_only = GetNonCstWorkFluidVelocity();
auto& vel_only = variable_holder.GetNonCstWorkFluidVelocity();
Wrappers::Petsc::MatMult(interpolator_holder.GetMonolithic2Velocity().GetInterpolationMatrix(),
GetDarcyOperator().GetVelocitySolutionParam().GetGlobalVector(),
......@@ -647,18 +652,18 @@ namespace HappyHeart
void VariationalFormulation<HyperelasticLawT>::SetSolidVelocityForDifferential(Vec petsc_solid_velocity)
{
auto& delta_displacement = GetNonCstDeltaDisplacement();
auto& work_vel_fluid = GetNonCstWorkFluidVelocity();
decltype(auto) variable_holder = this->GetNonCstVariableHolder();
auto& work_vel_fluid = variable_holder.GetNonCstWorkFluidVelocity();
delta_displacement.SetFromPetscVec(petsc_solid_velocity, __FILE__, __LINE__);
decltype(auto) variable_holder = this->GetVariableHolder();
{
variable_holder.template DevPrint<DevPhase::dH>(delta_displacement,
"initial_disp_delta");
}
Wrappers::Petsc::MatMult(variable_holder_parent::GetVariableHolder().GetSolidOnFluidMesh().GetMatrixSolidToFluid(),
delta_displacement,
......
Sources/ModelInstances/UnderDevelopment/Poromechanics/ImplicitStepFluid/InnerLoop/VariationalFormulationAccessors.hxx
View file @ 997dc1b1
......@@ -120,22 +120,7 @@ namespace HappyHeart
{
return const_cast<GlobalVector&>(GetWorkFluidMass());
}
template<class HyperelasticLawT>
inline const GlobalVector& VariationalFormulation<HyperelasticLawT>::GetWorkFluidVelocity() const noexcept
{
assert(!(!work_fluid_velocity_));
return *work_fluid_velocity_;
}
template<class HyperelasticLawT>
inline GlobalVector& VariationalFormulation<HyperelasticLawT>::GetNonCstWorkFluidVelocity() noexcept
{
return const_cast<GlobalVector&>(GetWorkFluidVelocity());
}
template<class HyperelasticLawT>
inline const GlobalVariationalOperatorNS::ScalarDivVectorial& VariationalFormulation<HyperelasticLawT>
......
Sources/ModelInstances/UnderDevelopment/Poromechanics/ImplicitStepFluid/InnerLoop/VariationalFormulationInit.hxx
View file @ 997dc1b1
......@@ -145,10 +145,6 @@ namespace HappyHeart
AllocateGlobalVector(god_of_dof, GetNonCstWorkFluidMass());
work_fluid_velocity_ =
std::make_unique<GlobalVector>(god_of_dof.GetNumberingSubset(EnumUnderlyingType(NumberingSubsetIndex::fluid_velocity)));
AllocateGlobalVector(god_of_dof, GetNonCstWorkFluidVelocity());
delta_fluid_mass_ = std::make_unique<GlobalVector>(GetWorkFluidMass());
......
Sources/ModelInstances/UnderDevelopment/Poromechanics/ImplicitStepFluid/VariationalFormulation.hxx
View file @ 997dc1b1
......@@ -78,6 +78,8 @@ namespace HappyHeart
double absolute_error = 1.e20;
double relative_error = 1.e20;
decltype(auto) interpolator_holder = this->GetInterpolatorHolder();
{
auto& velSHalfVhf = variable_holder.GetNonCstVelSHalfVhf();
......@@ -89,15 +91,17 @@ namespace HappyHeart
variable_holder.GetNonCstSolidOnFluidMesh().ComputeSolidVelocity(velSHalfVhf);
// \todo #820 Improve this!
// We get value on the whole mesh; we want to zero all values not on the Robin interface.
// Current stupid method is to interpolate back-and-forth; it's obviously not very efficient.
auto& velSHalfVhf_nullified_out_of_robin_interface
= variable_holder.GetNonCstVelSHalfVhfZeroedOutOfRobinInterface();
// We get value on the whole mesh; we want to zero all values not on the Robin interface.
Wrappers::Petsc::MatMult(interpolator_holder.GetVelocityNullifierOutOfRobinInterface(),
velSHalfVhf,
velSHalfVhf_nullified_out_of_robin_interface,
__FILE__, __LINE__);
variable_holder.template DevPrint<DevPhase::none>(velSHalfVhf,
"velSHalfVhfr_at_begin");
variable_holder.template DevPrint<DevPhase::none>(velSHalfVhf_nullified_out_of_robin_interface,
"velSHalfVhfr_non_robin_zero_at_begin");
}
......@@ -110,9 +114,6 @@ namespace HappyHeart
++newton_fp_index;
}
decltype(auto) interpolator_holder = this->GetInterpolatorHolder();
// Update fluidmass, velocity fluid and pressure.
Wrappers::Petsc::MatMult(interpolator_holder.GetMonolithic2Mass().GetInterpolationMatrix(),
inner_varf.GetCorrectedValues(),
......
Sources/ModelInstances/UnderDevelopment/Poromechanics/InterpolatorHolder.cpp
View file @ 997dc1b1
......@@ -134,10 +134,18 @@ namespace HappyHeart
std::make_unique<type>(fluid_felt_space,
fluid_god_of_dof.GetNumberingSubset(EnumUnderlyingType(NumberingSubsetIndex::fluid_velocity)),
robin_felt_space,
fluid_god_of_dof.GetNumberingSubset(EnumUnderlyingType(NumberingSubsetIndex::fluid_velocity)));
fluid_god_of_dof.GetNumberingSubset(EnumUnderlyingType(NumberingSubsetIndex::fluid_velocity_on_robin_interface)));
velocity_to_robin_interface_->Init();
velocity_robin_interface_to_full_mesh_ =
std::make_unique<type>(robin_felt_space,
fluid_god_of_dof.GetNumberingSubset(EnumUnderlyingType(NumberingSubsetIndex::fluid_velocity_on_robin_interface)),
fluid_felt_space,
fluid_god_of_dof.GetNumberingSubset(EnumUnderlyingType(NumberingSubsetIndex::fluid_velocity)));
velocity_robin_interface_to_full_mesh_->Init();
monolithic_2_robin_velocity_ =
std::make_unique<GlobalMatrix>(fluid_god_of_dof.GetNumberingSubset(EnumUnderlyingType(NumberingSubsetIndex::fluid_velocity)),
fluid_god_of_dof.GetNumberingSubset(EnumUnderlyingType(NumberingSubsetIndex::fluid_mass_velocity_pressure)));
......@@ -148,6 +156,22 @@ namespace HappyHeart
__FILE__, __LINE__);
}
{
decltype(auto) numbering_subset =
fluid_god_of_dof.GetNumberingSubset(EnumUnderlyingType(NumberingSubsetIndex::fluid_velocity));
velocity_nullifier_out_of_robin_interface_ = std::make_unique<GlobalMatrix>(numbering_subset,
numbering_subset);
Wrappers::Petsc::MatMatMult(GetVelocityRobinInterfaceToFullMesh().GetInterpolationMatrix(),
GetVelocityToRobinInterface().GetInterpolationMatrix(),
*velocity_nullifier_out_of_robin_interface_,
__FILE__, __LINE__);
}
}
......
Sources/ModelInstances/UnderDevelopment/Poromechanics/InterpolatorHolder.hpp
View file @ 997dc1b1
......@@ -113,6 +113,9 @@ namespace HappyHeart
//! Constant accessor to the matrix to go from (whole mesh, monolithic) to (robin interface, velocity).
const GlobalMatrix& GetMonolithic2RobinVelocity() const noexcept;
//! Constant accessor to the matrix used to zero all velocity terms that aren't on the Robin interface.
const GlobalMatrix& GetVelocityNullifierOutOfRobinInterface() const noexcept;
private:
/*!
......@@ -121,6 +124,7 @@ namespace HappyHeart
*/
GlobalMatrix& GetNonCstMonolithic2RobinVelocity() noexcept;
private:
//! Interpolator fluid mass on solid mesh -> fluid mass on fluid mesh.
......@@ -173,7 +177,8 @@ namespace HappyHeart
//! Matrix to go from (whole mesh, monolithic) to (robin interface, velocity).
GlobalMatrix::unique_ptr monolithic_2_robin_velocity_ = nullptr;
//! Matrix used to zero all velocity terms that aren't on the Robin interface.
GlobalMatrix::unique_ptr velocity_nullifier_out_of_robin_interface_ = nullptr;
///@}
......
Sources/ModelInstances/UnderDevelopment/Poromechanics/InterpolatorHolder.hxx
View file @ 997dc1b1
......@@ -131,6 +131,12 @@ namespace HappyHeart
}
inline const GlobalMatrix& InterpolatorHolder
::GetVelocityNullifierOutOfRobinInterface() const noexcept
{
assert(!(!velocity_nullifier_out_of_robin_interface_));
return *velocity_nullifier_out_of_robin_interface_;
}
} // namespace PoromechanicsNS
......
Sources/ModelInstances/UnderDevelopment/Poromechanics/VariableHolder.cpp
View file @ 997dc1b1
......@@ -115,6 +115,9 @@ namespace HappyHeart
fluid_god_of_dof.GetNumberingSubset(EnumUnderlyingType(NumberingSubsetIndex::fluid_velocity));
velSHalfVhf_ = std::make_unique<GlobalVector>(numbering_subset);
AllocateGlobalVector(fluid_god_of_dof, GetNonCstVelSHalfVhf());
work_fluid_velocity_ = std::make_unique<GlobalVector>(GetVelSHalfVhf());
velSHalfVhf_zeroed_out_of_robin_interface_ = std::make_unique<GlobalVector>(GetVelSHalfVhf());
}
}
......
Sources/ModelInstances/UnderDevelopment/Poromechanics/VariableHolder.hpp
View file @ 997dc1b1
......@@ -153,6 +153,15 @@ namespace HappyHeart
//! Non constant accessor to the odo #820 Rename this (Currently same as Bruno's Freefem)
GlobalVector& GetNonCstVelSHalfVhf() noexcept;
//! Constant accessor to the odo #820 Rename this (Currently same as Bruno's Freefem)
const GlobalVector& GetVelSHalfVhfZeroedOutOfRobinInterface() const noexcept;
//! Non constant accessor to the odo #820 Rename this (Currently same as Bruno's Freefem)
GlobalVector& GetNonCstVelSHalfVhfZeroedOutOfRobinInterface() noexcept;
//! Non constant accessor to the \#820 [DEV] Vector defined on fluid velocity numbering subset.
GlobalVector& GetNonCstWorkFluidVelocity() const noexcept;
//! Debug indexes...
unsigned int implicit_loop_index = 0u;
......@@ -170,6 +179,8 @@ namespace HappyHeart
// \todo #820 At the moment Bruno's name is kept!
GlobalVector::unique_ptr velSHalfVhf_ = nullptr;
GlobalVector::unique_ptr velSHalfVhf_zeroed_out_of_robin_interface_ = nullptr;
///@}
......@@ -216,7 +227,12 @@ namespace HappyHeart
//!
const TimeManager& time_manager_;
private:
//! \#820 [DEV] Vector defined on fluid velocity numbering subset.
mutable GlobalVector::unique_ptr work_fluid_velocity_ = nullptr;
};
......
Sources/ModelInstances/UnderDevelopment/Poromechanics/VariableHolder.hxx
View file @ 997dc1b1
......@@ -159,7 +159,27 @@ namespace HappyHeart
{
return const_cast<GlobalVector&>(GetVelSHalfVhf());
}
inline const GlobalVector& VariableHolder::GetVelSHalfVhfZeroedOutOfRobinInterface() const noexcept
{
assert(!(!velSHalfVhf_zeroed_out_of_robin_interface_));
return *velSHalfVhf_zeroed_out_of_robin_interface_;
}
inline GlobalVector& VariableHolder::GetNonCstVelSHalfVhfZeroedOutOfRobinInterface() noexcept
{
return const_cast<GlobalVector&>(GetVelSHalfVhfZeroedOutOfRobinInterface());
}
inline GlobalVector& VariableHolder::GetNonCstWorkFluidVelocity() const noexcept
{
assert(!(!work_fluid_velocity_));
return *work_fluid_velocity_;
}
} // namespace PoromechanicsNS
......
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