Skip to content

GitLab

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

#820 PerformDifferential was unduly in a manual Newton loop.

parent 8155078b
Hide whitespace changes
Inline Side-by-side
Sources/ModelInstances/UnderDevelopment/Poromechanics/ImplicitStepFluid/InnerLoop/VariationalFormulation.hxx
View file @ afc91187
......@@ -458,13 +458,10 @@ namespace HappyHeart
PETSC_VIEWER_ASCII_MATLAB);
}
// Wrappers::Petsc::AXPY(-1.,
// delta_sol,
// delta_displacement,
// __FILE__, __LINE__);
Wrappers::Petsc::AXPY(-1.,
delta_sol,
delta_displacement,
__FILE__, __LINE__);
{
const auto filename = parent::GetOutputDirectory(numbering_subset)
......
Sources/ModelInstances/UnderDevelopment/Poromechanics/ImplicitStepFluid/VariationalFormulation.hpp
View file @ afc91187
......@@ -131,9 +131,11 @@ namespace HappyHeart
* - Call Petsc's solver to solve the system.
* - Write the solutions.
*
* \param[in] is_differential If yes, run the differential of the variational formulation.
*/
void Perform(differential is_differential);
void Perform();
//! Perform differential step.
void PerformDifferential();
//! Constant accessor to the darcy vector after inner loop.
const GlobalVector& GetDarcyVector() const noexcept;
......
Sources/ModelInstances/UnderDevelopment/Poromechanics/ImplicitStepFluid/VariationalFormulation.hxx
View file @ afc91187
......@@ -60,7 +60,7 @@ namespace HappyHeart
template<class HyperelasticLawT>
void VariationalFormulation<HyperelasticLawT>::Perform(differential is_differential)
void VariationalFormulation<HyperelasticLawT>::Perform()
{
double current_norm = 1.e20;
......@@ -71,89 +71,81 @@ namespace HappyHeart
auto counter = 0u;
// Newton loop (\todo #820 To replace with Petsc Newton...)
do
current_norm = inner_varf.Perform(counter++);
// Update fluidmass, velocity fluid and pressure.
Wrappers::Petsc::MatMult(inner_varf.GetMonolithic2Mass().GetInterpolationMatrix(),
inner_varf.GetCorrectedValues(),
variable_holder_parent::GetNonCstVariableHolder().GetNonCstFluidMassVector(),
__FILE__, __LINE__);
// Perform one more time to update T21, T33 and T33 without Dirichlet.
{
switch(is_differential)
{
case differential::no:
{
current_norm = inner_varf.Perform(counter++);
break;
}
case differential::yes:
{
current_norm = inner_varf.PerformDifferential(counter++,
GetMonolithicMatrixAfterInnerLoop(),
GetMonolithicMatrixAfterInnerLoopWithoutDirichlet());
break;
}
}
decltype(auto) updated_matrix = inner_varf.ComputeMatrixWithUpdatedT21();
std::cout << "Norm after inner loop = " << current_norm << std::endl;
} while (current_norm > 1.e-12 && counter < 10u); // \todo #820 Replace by proper stop condition from input file!
// \todo #820 - #530 Use Swap when available!
GetNonCstMonolithicMatrixAfterInnerLoop().Copy(updated_matrix, __FILE__, __LINE__);
updated_matrix.View(parent::MpiHappyHeart(),
parent::GetOutputDirectory(numbering_subset)
+ "/"
+ DifferentialPreffix(differential::no)
+ "matrix_after_FP_mixt_with_dirichlet_time_"
+ std::to_string(parent::GetTimeManager().NtimeModified()) + ".m",
__FILE__, __LINE__,
PETSC_VIEWER_ASCII_MATLAB);
}
if (is_differential == differential::no)
{
// Update fluidmass, velocity fluid and pressure.
Wrappers::Petsc::MatMult(inner_varf.GetMonolithic2Mass().GetInterpolationMatrix(),
inner_varf.GetCorrectedValues(),
variable_holder_parent::GetNonCstVariableHolder().GetNonCstFluidMassVector(),
__FILE__, __LINE__);
// \todo #820 - #530 Use Swap when available!
// Perform one more time to update T21, T33 and T33 without Dirichlet.
{
decltype(auto) updated_matrix = inner_varf.ComputeMatrixWithUpdatedT21();
// \todo #820 - #530 Use Swap when available!
GetNonCstMonolithicMatrixAfterInnerLoop().Copy(updated_matrix, __FILE__, __LINE__);
updated_matrix.View(parent::MpiHappyHeart(),
parent::GetOutputDirectory(numbering_subset)
+ "/"
+ DifferentialPreffix(is_differential)
+ "matrix_after_FP_mixt_with_dirichlet_time_"
+ std::to_string(parent::GetTimeManager().NtimeModified()) + ".m",
__FILE__, __LINE__,
PETSC_VIEWER_ASCII_MATLAB);
}
decltype(auto) updated_matrix = inner_varf.ComputeMatrixWithUpdatedT33WithoutDirichlet();
GetNonCstMonolithicMatrixAfterInnerLoopWithoutDirichlet().Copy(updated_matrix, __FILE__, __LINE__);
{
// \todo #820 - #530 Use Swap when available!
decltype(auto) updated_matrix = inner_varf.ComputeMatrixWithUpdatedT33WithoutDirichlet();
GetNonCstMonolithicMatrixAfterInnerLoopWithoutDirichlet().Copy(updated_matrix, __FILE__, __LINE__);
// \todo #820 Not exactly equal to Freefem; investigate if one term is missing!
updated_matrix.View(parent::MpiHappyHeart(),
parent::GetOutputDirectory(numbering_subset)
+ "/"
+ DifferentialPreffix(is_differential)
+ "matrix_after_FP_mixt_without_dirichlet_time_"
+ std::to_string(parent::GetTimeManager().NtimeModified()) + ".m",
__FILE__, __LINE__,
PETSC_VIEWER_ASCII_MATLAB);
}
// \todo #820 Not exactly equal to Freefem; investigate if one term is missing!
updated_matrix.View(parent::MpiHappyHeart(),
parent::GetOutputDirectory(numbering_subset)
+ "/"
+ DifferentialPreffix(differential::no)
+ "matrix_after_FP_mixt_without_dirichlet_time_"
+ std::to_string(parent::GetTimeManager().NtimeModified()) + ".m",
__FILE__, __LINE__,
PETSC_VIEWER_ASCII_MATLAB);
}
{
decltype(auto) darcy_vector = GetNonCstDarcyVector();
{
decltype(auto) darcy_vector = GetNonCstDarcyVector();
inner_varf.ComputeDarcy(-1., IsFullDarcy::no, darcy_vector);
darcy_vector.View(parent::MpiHappyHeart(),
parent::GetOutputDirectory(numbering_subset)
+ "/"
+ DifferentialPreffix(is_differential)
+ "darcy_vector_after_FP_mixt_time_"
+ std::to_string(parent::GetTimeManager().NtimeModified()) + ".m",
__FILE__, __LINE__,
PETSC_VIEWER_ASCII_MATLAB);
}
inner_varf.ComputeDarcy(-1., IsFullDarcy::no, darcy_vector);
darcy_vector.View(parent::MpiHappyHeart(),
parent::GetOutputDirectory(numbering_subset)
+ "/"
+ DifferentialPreffix(differential::no)
+ "darcy_vector_after_FP_mixt_time_"
+ std::to_string(parent::GetTimeManager().NtimeModified()) + ".m",
__FILE__, __LINE__,
PETSC_VIEWER_ASCII_MATLAB);
}
}
template<class HyperelasticLawT>
void VariationalFormulation<HyperelasticLawT>::PerformDifferential()
{
double current_norm = 1.e20;
auto& inner_varf = GetNonCstInnerLoopVarf();
auto counter = 0u;
current_norm = inner_varf.PerformDifferential(counter++,
GetMonolithicMatrixAfterInnerLoop(),
GetMonolithicMatrixAfterInnerLoopWithoutDirichlet());
std::cout << "Norm after inner loop = " << current_norm << std::endl;
}
template<class HyperelasticLawT>
void VariationalFormulation<HyperelasticLawT>::AllocateMatricesAndVectors()
{
......
Sources/ModelInstances/UnderDevelopment/Poromechanics/Model.hxx
View file @ afc91187
......@@ -48,7 +48,7 @@ namespace HappyHeart
}
auto& implicit_step_fluid_varf = GetNonCstImplicitStepFluidVarf();
implicit_step_fluid_varf.Perform(differential::no);
implicit_step_fluid_varf.Perform();
UpdateSolidVector();
......@@ -247,7 +247,7 @@ namespace HappyHeart
auto& implicit_step_varf = GetNonCstImplicitStepFluidVarf();
implicit_step_varf.SetSolidVelocityForDifferential(petsc_input_vector);
implicit_step_varf.Perform(differential::yes);
implicit_step_varf.PerformDifferential();
// We can't take directly the internal pointer of a wrapped vector which pattern would be similar to input_displ_vector
......
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