Skip to content

GitLab

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

#820 Investigating the solve issue in second implicit loop iteration.

parent 4a82bfc9
Hide whitespace changes
Inline Side-by-side
Sources/ModelInstances/Hyperelasticity/Private/VectorsAndMatricesCrtp.hpp
View file @ c41fa25c
......@@ -84,12 +84,29 @@ namespace HappyHeart
/*!
* \brief Update the vector that contains the values seeked at the moment the residual is evaluated.
* \brief Update the vector that contains the values sought at the moment the residual is evaluated.
*
* \param[in] evaluation_state Evaluation state given by the Petsc's snes function.
*/
void UpdateEvaluationState(Vec evaluation_state);
/*!
* \brief Set EvaluationState value from a given global vector.
*
* This should not be used during a Newton loop except in the initialization: if you need a non zero
* starting point you may use it through something like:
* \code
* if (this->GetSnes().GetSnesIteration(__FILE__, __LINE__) > 0)
* vm.UpdateEvaluationState(a_evaluation_state); // Vec given by Petsc
* else
* vm.SetInitialEvaluationState(my_evaluation_state); // GlobalVector given manually.
* \endcode
*
* \param[in] evaluation_state Evaluation state
*/
void SetInitialEvaluationState(const GlobalVector& evaluation_state,
const char* invoking_file, int invoking_line);
/// \name Accessors to the global vectors and matrices managed by the class.
///@{
......
Sources/ModelInstances/Hyperelasticity/Private/VectorsAndMatricesCrtp.hxx
View file @ c41fa25c
......@@ -53,6 +53,16 @@ namespace HappyHeart
}
template<class DerivedT>
void VectorsAndMatricesCrtp<DerivedT>
::SetInitialEvaluationState(const GlobalVector& value,
const char* invoking_file, int invoking_line)
{
auto& evaluation_state = GetNonCstEvaluationState();
evaluation_state.Copy(value, invoking_file, invoking_line);
}
template<class DerivedT>
inline const GlobalVector& VectorsAndMatricesCrtp<DerivedT>::GetForce() const noexcept
{
......
Sources/ModelInstances/UnderDevelopment/Poromechanics/Model.hpp
View file @ c41fa25c
......@@ -224,14 +224,7 @@ namespace HappyHeart
const implicit_step_fluid_varf_type& GetImplicitStepFluidVarf() const noexcept;
//! Non constant accessor to the variational formulation for fluid implicit step
implicit_step_fluid_varf_type& GetNonCstImplicitStepFluidVarf() noexcept;
//! Access to fluid mass parameter.
const ParameterAtDof<ParameterNS::Type::scalar>::type& GetFluidMassParameter() const noexcept;
//! Parameter that encapsulates fluid_mass_vector from previous time iteration, expressed on the solid mesh.
const ParameterAtDof<ParameterNS::Type::scalar, ParameterNS::TimeDependencyNS::None>::type&
GetPreviousTimeFluidMassParameter() const noexcept;
implicit_step_fluid_varf_type& GetNonCstImplicitStepFluidVarf() noexcept;
//! Constant accessor to the computed solid displacement.
const GlobalVector& GetSolidDisplacement() const noexcept;
......@@ -428,13 +421,7 @@ namespace HappyHeart
//! Object that includes the interpolators.
InterpolatorHolder::const_unique_ptr interpolator_holder_ = nullptr;
//! Parameter that encapsulates fluid_mass_vector, expressed on the solid mesh.
ParameterAtDof<ParameterNS::Type::scalar>::type::unique_ptr fluid_mass_param_ = nullptr;
//! Parameter that encapsulates fluid_mass_vector from previous time iteration, expressed on the solid mesh.
ParameterAtDof<ParameterNS::Type::scalar, ParameterNS::TimeDependencyNS::None>::type::unique_ptr
prev_time_it_fluid_mass_param_ = nullptr;
/*!
* \brief Initial porosity value as read in the input parameter file.
*
......
Sources/ModelInstances/UnderDevelopment/Poromechanics/Model.hxx
View file @ c41fa25c
......@@ -39,11 +39,8 @@ namespace HappyHeart
auto& iteration_index = variable_holder.implicit_loop_index;
iteration_index = 0u;
auto& solid_displacement = variable_holder.GetNonCstSolidDisplacement();
while (absolute_error > 1.e-10 && relative_error > 1.e-6 && iteration_index < 6u) // \todo #820 Should not remain hardcoded!
// \todo #820 iteration index sup value is obviously temporary for dev purposes...
{
......
Sources/ModelInstances/UnderDevelopment/Poromechanics/ModelAccessors.hxx
View file @ c41fa25c
......@@ -111,26 +111,6 @@ namespace HappyHeart
}
template<class SolidVariationalFormulationPolicyT>
inline const ParameterAtDof<ParameterNS::Type::scalar>::type& Model<SolidVariationalFormulationPolicyT>
::GetFluidMassParameter() const noexcept
{
assert(!(!fluid_mass_param_));
return *fluid_mass_param_;
}
template<class SolidVariationalFormulationPolicyT>
inline const ParameterAtDof<ParameterNS::Type::scalar, ParameterNS::TimeDependencyNS::None>::type&
Model<SolidVariationalFormulationPolicyT>
::GetPreviousTimeFluidMassParameter() const noexcept
{
assert(!(!prev_time_it_fluid_mass_param_));
return *prev_time_it_fluid_mass_param_;
}
template<class SolidVariationalFormulationPolicyT>
inline const GlobalVector& Model<SolidVariationalFormulationPolicyT>::GetSolidDisplacement() const noexcept
{
......
Sources/ModelInstances/UnderDevelopment/Poromechanics/ModelInitialize.hxx
View file @ c41fa25c
......@@ -114,7 +114,6 @@ namespace HappyHeart
decltype(auto) solid_god_of_dof = parent::GetGodOfDof(EnumUnderlyingType(MeshIndex::solid));
decltype(auto) fluid_mesh = fluid_god_of_dof.GetGeometricMeshRegion();
decltype(auto) solid_mesh = solid_god_of_dof.GetGeometricMeshRegion();
decltype(auto) unknown_manager = UnknownManager::GetInstance();
{
decltype(auto) numbering_subset =
......@@ -134,29 +133,7 @@ namespace HappyHeart
{
const auto& felt_space = solid_god_of_dof.GetFEltSpace(EnumUnderlyingType(FEltSpaceIndex::solid));
using scalar_at_dof_type = ParameterAtDof<ParameterNS::Type::scalar>::type;
const auto& fluid_mass_unknown = unknown_manager.GetUnknown(EnumUnderlyingType(UnknownIndex::fluid_mass));
fluid_mass_param_ = std::make_unique<scalar_at_dof_type>("Fluid mass",
solid_mesh,
felt_space,
fluid_mass_unknown,
GetVariableHolder().GetFluidMassVector());
decltype(auto) inlet_border_felt_space
= solid_god_of_dof.GetFEltSpace(EnumUnderlyingType(FEltSpaceIndex::solid_inlet_border));
prev_time_it_fluid_mass_param_ = std::make_unique<scalar_at_dof_type>("Fluid mass previous time iteration on inlet border",
solid_mesh,
inlet_border_felt_space,
fluid_mass_unknown,
GetVariableHolder().GetFluidMassVectorPreviousIteration());
}
InitPorosity(input_parameter_data);
......@@ -188,7 +165,7 @@ namespace HappyHeart
const auto& time_manager = parent::GetTimeManager();
const auto& mpi = parent::MpiHappyHeart();
decltype(auto) fluid_density = GetFluidDensity();
decltype(auto) variable_holder = this->GetVariableHolder();
decltype(auto) solid_god_of_dof = parent::GetGodOfDof(EnumUnderlyingType(MeshIndex::solid));
......@@ -202,7 +179,7 @@ namespace HappyHeart
numbering_subset,
time_manager,
solid_god_of_dof,
GetFluidMassParameter(),
variable_holder.GetFluidMassParameter(),
fluid_density,
GetSolidDisplacement());
......@@ -240,6 +217,7 @@ namespace HappyHeart
decltype(auto) bc_manager = ::HappyHeart::DirichletBoundaryConditionManager::GetInstance();
decltype(auto) fluid_density = GetFluidDensity();
decltype(auto) fluid_viscosity = GetFluidViscosity();
decltype(auto) variable_holder = this->GetVariableHolder();
{
decltype(auto) numbering_subset =
......@@ -275,7 +253,7 @@ namespace HappyHeart
fluid_density,
GetPenalizationPorosity(),
GetInitialPorosityValue(),
GetFluidMassParameter(),
variable_holder.GetFluidMassParameter(),
GetBulkSolid());
}
......@@ -392,13 +370,14 @@ namespace HappyHeart
decltype(auto) felt_space = god_of_dof.GetFEltSpace(EnumUnderlyingType(FEltSpaceIndex::solid_inlet_border));
decltype(auto) solid_displacement =
UnknownManager::GetInstance().GetUnknown(EnumUnderlyingType(UnknownIndex::solid_displacement));
decltype(auto) variable_holder = this->GetVariableHolder();
load_on_solid_operator_ =
std::make_unique<GlobalVariationalOperatorNS::LoadOnSolid>(felt_space,
solid_displacement,
mesh.GetDimension(),
GetInletPressure(),
GetPreviousTimeFluidMassParameter(),
variable_holder.GetPreviousTimeFluidMassParameter(),
GetFluidDensity(),
DoComputeProcessorWiseLocal2Global::yes // \todo #820 no is probably ok.
);
......
Sources/ModelInstances/UnderDevelopment/Poromechanics/SolidVariationalFormulationPolicy/Hyperelasticity/VariationalFormulation.hxx
View file @ c41fa25c
......@@ -234,12 +234,18 @@ namespace HappyHeart
const auto& displacement_previous_time_iteration =
vectors_and_matrices.GetDisplacementPreviousTimeIteration();
variable_holder.template DevPrint<DevPhase::solid_newton>(vectors_and_matrices.GetMidpointPosition(),
"dispSInterr_before_solid_solve");
variable_holder.template DevPrint<DevPhase::solid_newton>(vectors_and_matrices.GetEvaluationState(),
"dispSr_before_solid_solve");
variable_holder.template DevPrint<DevPhase::solid_newton>(displacement_previous_time_iteration,
"dispSr_old_before_solid_solve");
variable_holder.template DevPrint<DevPhase::solid_newton>(variable_holder.GetFluidMassParameter().GetGlobalVector(),
"fluidmass_before_solid_solve");
const auto& numbering_subset = GetNumberingSubset();
......@@ -248,7 +254,8 @@ namespace HappyHeart
HyperelasticityNS::Private::ComputeDynamicTimeSchemeResidualContribution(vectors_and_matrices, rhs);
variable_holder.template DevPrint<DevPhase::solid_newton>(rhs,
"solid_solve_residual_0");
const auto& velocity_previous_time_iteration =
vectors_and_matrices.GetVelocityPreviousTimeIteration();
......@@ -268,13 +275,23 @@ namespace HappyHeart
auto& dynamic_contribution_to_rhs = this->GetNonCstHelperGlobalVector<1>();
Wrappers::Petsc::MatMult(vectors_and_matrices.GetMassPerSquareTime(), helper_vector,
Wrappers::Petsc::MatMult(vectors_and_matrices.GetMassPerSquareTime(),
helper_vector,
dynamic_contribution_to_rhs, __FILE__, __LINE__);
Wrappers::Petsc::AXPY(1., dynamic_contribution_to_rhs, rhs, __FILE__, __LINE__);
Wrappers::Petsc::AXPY(1.,
dynamic_contribution_to_rhs,
rhs,
__FILE__, __LINE__);
variable_holder.template DevPrint<DevPhase::solid_newton>(rhs,
"solid_solve_residual_1");
dof_source_parent::AddToRhs(rhs);
variable_holder.template DevPrint<DevPhase::solid_newton>(rhs,
"solid_solve_residual_2");
parent::template ApplyEssentialBoundaryCondition
<
OperatorNS::Nature::linear,
......@@ -335,7 +352,17 @@ namespace HappyHeart
{
auto& vm = this->GetNonCstVectorsAndMatrices();
vm.UpdateEvaluationState(a_evaluation_state);
if (this->GetSnes().GetSnesIteration(__FILE__, __LINE__) > 0)
{
std::cout << "USUAL CASE" << std::endl;
vm.UpdateEvaluationState(a_evaluation_state);
}
else
{
std::cout << "Init CASE" << std::endl;
decltype(auto) variable_holder = this->GetNonCstVariableHolder();
vm.SetInitialEvaluationState(variable_holder.GetSolidDisplacement(), __FILE__, __LINE__);
}
{
// Compute some matrices and vectors through arithmetical operations (e.g. vector for midpoint position).
......
Sources/ModelInstances/UnderDevelopment/Poromechanics/VariableHolder.cpp
View file @ c41fa25c
......@@ -35,6 +35,7 @@ namespace HappyHeart
decltype(auto) solid_god_of_dof = god_of_dof_manager.GetGodOfDof(EnumUnderlyingType(MeshIndex::solid));
decltype(auto) fluid_mesh = fluid_god_of_dof.GetGeometricMeshRegion();
decltype(auto) solid_mesh = solid_god_of_dof.GetGeometricMeshRegion();
decltype(auto) unknown_manager = UnknownManager::GetInstance();
{
namespace ipl = Utilities::InputParameterListNS;
......@@ -135,6 +136,31 @@ namespace HappyHeart
}
{
const auto& felt_space = solid_god_of_dof.GetFEltSpace(EnumUnderlyingType(FEltSpaceIndex::solid));
using scalar_at_dof_type = ParameterAtDof<ParameterNS::Type::scalar>::type;
const auto& fluid_mass_unknown = unknown_manager.GetUnknown(EnumUnderlyingType(UnknownIndex::fluid_mass));
fluid_mass_param_ = std::make_unique<scalar_at_dof_type>("Fluid mass",
solid_mesh,
felt_space,
fluid_mass_unknown,
GetFluidMassVector());
decltype(auto) inlet_border_felt_space
= solid_god_of_dof.GetFEltSpace(EnumUnderlyingType(FEltSpaceIndex::solid_inlet_border));
prev_time_it_fluid_mass_param_ = std::make_unique<scalar_at_dof_type>("Fluid mass previous time iteration on inlet border",
solid_mesh,
inlet_border_felt_space,
fluid_mass_unknown,
GetFluidMassVectorPreviousIteration());
}
}
......
Sources/ModelInstances/UnderDevelopment/Poromechanics/VariableHolder.hpp
View file @ c41fa25c
......@@ -21,6 +21,7 @@
# include "Parameters/Parameter.hpp"
# include "Parameters/Compound/Solid/Solid.hpp"
# include "Parameters/ParameterAtDof.hpp"
# include "ModelInstances/UnderDevelopment/Poromechanics/Private/SolidOnFluidMesh.hpp"
......@@ -178,6 +179,15 @@ namespace HappyHeart
unsigned int dh_index = 0u;
unsigned int solid_newton_index = 0u;
//! Access to fluid mass parameter.
const ParameterAtDof<ParameterNS::Type::scalar>::type& GetFluidMassParameter() const noexcept;
//! Parameter that encapsulates fluid_mass_vector from previous time iteration on the inlet border,
//! expressed on the solid mesh.
const ParameterAtDof<ParameterNS::Type::scalar, ParameterNS::TimeDependencyNS::None>::type&
GetPreviousTimeFluidMassParameter() const noexcept;
private:
......@@ -236,6 +246,15 @@ namespace HappyHeart
//!
const TimeManager& time_manager_;
//! Parameter that encapsulates fluid_mass_vector, expressed on the solid mesh.
ParameterAtDof<ParameterNS::Type::scalar>::type::unique_ptr fluid_mass_param_ = nullptr;
//! Parameter that encapsulates fluid_mass_vector from previous time iteration, expressed on the solid mesh.
ParameterAtDof<ParameterNS::Type::scalar, ParameterNS::TimeDependencyNS::None>::type::unique_ptr
prev_time_it_fluid_mass_param_ = nullptr;
private:
......
Sources/ModelInstances/UnderDevelopment/Poromechanics/VariableHolder.hxx
View file @ c41fa25c
......@@ -195,6 +195,25 @@ namespace HappyHeart
}
inline const ParameterAtDof<ParameterNS::Type::scalar>::type& VariableHolder
::GetFluidMassParameter() const noexcept
{
assert(!(!fluid_mass_param_));
return *fluid_mass_param_;
}
inline const ParameterAtDof<ParameterNS::Type::scalar, ParameterNS::TimeDependencyNS::None>::type&
VariableHolder
::GetPreviousTimeFluidMassParameter() const noexcept
{
assert(!(!prev_time_it_fluid_mass_param_));
return *prev_time_it_fluid_mass_param_;
}
} // 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