Skip to content

GitLab

Commit 3cd44f15 authored by GILLES Sebastien's avatar GILLES Sebastien
Browse files

#723 FSI/Newton: rewrite a bit so that displacement used to compute new... 

#723 FSI/Newton: rewrite a bit so that displacement used to compute new velocity applied on Dirichket is an argument.
parent 53cfd748
Hide whitespace changes
Inline Side-by-side
Sources/ModelInstances/FSI_EI/Newton/Model.hpp
View file @ 3cd44f15
......@@ -156,8 +156,19 @@ namespace HappyHeart
//! Implicit fluid step.
void ImplicitFluidStep(differential is_differential);
//! Update fluid Dirichlet condition.
void UpdateFluidDirichletCondition(differential is_differential);
/*!
* \brief Update fluid Dirichlet condition.
*
* \internal Must be called immediately after solid_velocity_on_interface_ has been updated; it is not
* performed within the method because differential and non differential cases use different values
* for this velocity.
*/
void UpdateFluidDirichletCondition();
void UpdateSolidVelocityOnInterfaceAfterFluidImplicitStep();
void UpdateSolidVelocityOnInterfaceAfterDiffFluidImplicitStep(const GlobalVector& solid_displacement);
//! Update solid displacements before SoF.
void UpdateSolidDisplacementBeforeSoF();
......@@ -232,6 +243,12 @@ namespace HappyHeart
//! Solid displacement before SoF is applied (in Aitken loop). See Bruno's note for more details about Aitken.
GlobalVector& GetNonCstSolidDisplacementBeforeSoFOnInterface() noexcept;
//! Solid velocity on interface.
const GlobalVector& GetSolidVelocityOnInterface() const noexcept;
//! 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;
......@@ -259,6 +276,8 @@ namespace HappyHeart
GlobalVector::unique_ptr solid_displacement_after_SoF_ = nullptr;
GlobalVector::unique_ptr solid_displacement_after_dH_ = nullptr; // \todo #723 Add accessors if kept
GlobalVector::unique_ptr solid_displacement_evolution_in_SoF_ = nullptr;
GlobalVector::unique_ptr solid_displacement_evolution_in_SoF_on_interface_ = nullptr;
......
Sources/ModelInstances/FSI_EI/Newton/Model.hxx
View file @ 3cd44f15
......@@ -47,7 +47,6 @@ namespace HappyHeart
auto&& bc_list =
{ bc_manager.GetDirichletBoundaryConditionPtr("fluid_first"),
bc_manager.GetDirichletBoundaryConditionPtr("fluid_second") };
explicit_step_variational_formulation_ =
std::make_unique<ExplicitStepVariationalFormulation>(mpi,
......@@ -60,7 +59,6 @@ namespace HappyHeart
formulation.Init(input_parameter_data);
sum_fluid_residual_ = std::make_unique<GlobalVector>(formulation.GetSystemRhs(numbering_subset));
}
{
......@@ -90,7 +88,6 @@ namespace HappyHeart
}
{
const auto& solid_displacement = unknown_manager.GetUnknown(EnumUnderlyingType(UnknownIndex::solid_displacement));
const auto& main_felt_space = god_of_dof.GetFEltSpace(EnumUnderlyingType(FEltSpaceIndex::solid));
......@@ -137,6 +134,9 @@ namespace HappyHeart
solid_displacement_after_SoF_ =
std::make_unique<GlobalVector>(solid_displacement_before_SoF);
solid_displacement_after_dH_ =
std::make_unique<GlobalVector>(solid_displacement_before_SoF);
solid_displacement_previous_time_iteration_ =
std::make_unique<GlobalVector>(solid_displacement_before_SoF);
......
Sources/ModelInstances/FSI_EI/Newton/ModelAccessors.hxx
View file @ 3cd44f15
......@@ -167,7 +167,23 @@ namespace HappyHeart
return *differential_solid_solver_;
}
template<class SolidVariationalFormulationPolicyT>
const GlobalVector& Model<SolidVariationalFormulationPolicyT>::GetSolidVelocityOnInterface() const noexcept
{
assert(!(!solid_velocity_on_interface_));
return *solid_velocity_on_interface_;
}
template<class SolidVariationalFormulationPolicyT>
GlobalVector& Model<SolidVariationalFormulationPolicyT>::GetNonCstSolidVelocityOnInterface() noexcept
{
return const_cast<GlobalVector&>(GetSolidVelocityOnInterface());
}
} // namespace FSI_EINS
......
Sources/ModelInstances/FSI_EI/Newton/ModelForward.hxx
View file @ 3cd44f15
......@@ -36,7 +36,8 @@ namespace HappyHeart
UpdateSolidDisplacementBeforeSoF();
UpdateFluidDirichletCondition(differential::no);
UpdateSolidVelocityOnInterfaceAfterFluidImplicitStep();
UpdateFluidDirichletCondition();
ImplicitFluidStep(differential::no);
ComputeSolidResidual();
......@@ -45,9 +46,15 @@ namespace HappyHeart
UpdateSolidDisplacementAfterSoF();
UpdateFluidDirichletCondition(differential::yes);
auto& solid_displacement_after_dH = *solid_displacement_after_dH_;
solid_displacement_after_dH_->Copy(GetSolidDisplacementAfterSoF(), __FILE__, __LINE__);
UpdateSolidVelocityOnInterfaceAfterDiffFluidImplicitStep(GetSolidDisplacementEvolutionInSoFOnInterface());
UpdateFluidDirichletCondition();
ImplicitFluidStep(differential::yes);
auto& fluid_varf = GetNonCstImplicitStepFluidVariationalFormulation();
{
......@@ -116,44 +123,19 @@ namespace HappyHeart
<< " on " << solution.GetProgramWiseSize(__FILE__, __LINE__) << " elements." << std::endl; // opposite of the Freefem output
}
Wrappers::Petsc::AXPY(-1.,
solution,
solid_displacement_after_dH,
__FILE__, __LINE__);
{
const PetscReal evaluation_state_min = solid_displacement_after_dH.Min(__FILE__, __LINE__).second;
const PetscReal evaluation_state_max = solid_displacement_after_dH.Max(__FILE__, __LINE__).second;
std::cout << "corrected displacement EXTREMA ARE " << evaluation_state_min << " and " << evaluation_state_max
<< " on " << solid_displacement_after_dH.GetProgramWiseSize(__FILE__, __LINE__) << " elements." << std::endl; // opposite of the Freefem output
}
//
// const auto& sol = fluid_varf.GetSystemSolution(fluid_varf.GetNumberingSubset());
// sol.View(mpi, __FILE__, __LINE__);
// do
// {
// const std::string aitken_output_dir = output_dir + "/aitken_" + std::to_string(aitken_index) + "/";
// if (mpi.IsRootProcessor())
// Folder::Create(aitken_output_dir);
//
// mpi.Barrier();
//
// (aitken_output_dir);
//
// UpdateFluidDirichletCondition(aitken_output_dir);
//
// ImplicitFluidStep(aitken_output_dir, aitken_index);
//
// ComputeSolidResidual(aitken_output_dir);
//
// ImplicitSolidStep(aitken_output_dir, aitken_index);
//
// UpdateSolidDisplacementAfterSoF(aitken_output_dir);
//
// const double aitken_coefficient = ComputeAitkenCoefficient(aitken_index, omega_0);
//
// UpdateSolidDisplacementEndOfAitken(aitken_coefficient, aitken_output_dir);
//
// aitken_error = ComputeAitkenError();
//
// } while (aitken_error > aitken_tolerance && aitken_index++ < 100u);
}
......@@ -192,54 +174,58 @@ namespace HappyHeart
template<class SolidVariationalFormulationPolicyT>
void Model<SolidVariationalFormulationPolicyT>
::UpdateFluidDirichletCondition(differential is_differential)
::UpdateFluidDirichletCondition()
{
// Update the value of the Dirichlet condition.
// Compute first the solid velocity from the previous aitken iteration.
const double inv_time_step = 1. / this->GetTimeManager().GetTimeStep();
const auto& solid_velocity_on_interface = GetSolidVelocityOnInterface();
const auto& p1_2_p2_interpolation_matrix = solid_to_fluid_velocity_on_fsi_interpolator_->GetInterpolationMatrix();
auto& solid_velocity_on_interface = *solid_velocity_on_interface_; // \todo #608 Proper accessor!
auto& p2_solid_velocity_on_interface = *p2_solid_velocity_on_interface_; // \todo #608 Proper accessor!
auto& implicit_fluid_boundary_condition =
Private::DirichletBoundaryConditionManager::GetInstance().GetNonCstDirichletBoundaryCondition(EnumUnderlyingType(BoundaryConditionIndex::fluid_first));
switch(is_differential)
{
case differential::no:
{
const auto& solid_displacement_previous_time_iteration_on_interface = *solid_displacement_previous_time_iteration_on_interface_;
solid_velocity_on_interface.Copy(GetSolidDisplacementBeforeSoFOnInterface(),
__FILE__, __LINE__);
Wrappers::Petsc::AXPY(-1.,
solid_displacement_previous_time_iteration_on_interface,
solid_velocity_on_interface,
__FILE__, __LINE__);
break;
}
case differential::yes:
{
solid_velocity_on_interface.Copy(GetSolidDisplacementEvolutionInSoFOnInterface(),
__FILE__, __LINE__);
break;
}
}
solid_velocity_on_interface.Scale(inv_time_step, __FILE__, __LINE__);
// Then interpolate it to P2 space.
Wrappers::Petsc::MatMult(p1_2_p2_interpolation_matrix,
Wrappers::Petsc::MatMult(solid_to_fluid_velocity_on_fsi_interpolator_->GetInterpolationMatrix(),
solid_velocity_on_interface,
p2_solid_velocity_on_interface,
__FILE__, __LINE__);
// And update Dirichlet boundary condition.
auto& implicit_fluid_boundary_condition =
Private::DirichletBoundaryConditionManager::GetInstance().GetNonCstDirichletBoundaryCondition(EnumUnderlyingType(BoundaryConditionIndex::fluid_first));
implicit_fluid_boundary_condition.UpdateValues(p2_solid_velocity_on_interface);
}
template<class SolidVariationalFormulationPolicyT>
void Model<SolidVariationalFormulationPolicyT>
::UpdateSolidVelocityOnInterfaceAfterFluidImplicitStep()
{
const auto& solid_displacement_previous_time_iteration_on_interface = *solid_displacement_previous_time_iteration_on_interface_;
auto& solid_velocity_on_interface = GetNonCstSolidVelocityOnInterface();
solid_velocity_on_interface.Copy(GetSolidDisplacementBeforeSoFOnInterface(),
__FILE__, __LINE__);
Wrappers::Petsc::AXPY(-1.,
solid_displacement_previous_time_iteration_on_interface,
solid_velocity_on_interface,
__FILE__, __LINE__);
const double inv_time_step = 1. / this->GetTimeManager().GetTimeStep();
solid_velocity_on_interface.Scale(inv_time_step, __FILE__, __LINE__);
}
template<class SolidVariationalFormulationPolicyT>
void Model<SolidVariationalFormulationPolicyT>
::UpdateSolidVelocityOnInterfaceAfterDiffFluidImplicitStep(const GlobalVector& solid_displacement)
{
auto& solid_velocity_on_interface = GetNonCstSolidVelocityOnInterface();
solid_velocity_on_interface.Copy(solid_displacement,
__FILE__, __LINE__);
const double inv_time_step = 1. / this->GetTimeManager().GetTimeStep();
solid_velocity_on_interface.Scale(inv_time_step, __FILE__, __LINE__);
}
......
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