Commit 8155078b authored by GILLES Sebastien's avatar GILLES Sebastien
Browse files

#820 Refine indexes used to tag iteration by splitting explicitly...

#820 Refine indexes used to tag iteration by splitting explicitly differential/non differential case (useful for debug comparison with Freefem).
parent eb9444c0
......@@ -24,7 +24,7 @@ namespace HappyHeart
//! Enum class to choose whether it is run in differential mode or not.
enum class differential { no, yes };
enum class differential { no = 0, yes };
//! Returns preffix "differential_" if differential is yes.
......
......@@ -483,6 +483,7 @@ namespace HappyHeart
GlobalMatrix& GetNonCstT21OnFluidMesh() noexcept;
//! Index of the current inner loop (for outputs).
template<differential is_differential>
unsigned int GetIndexInnerLoop() const noexcept;
//! Constant accessor to the parameter that encapsulates velocity part of the current solution.
......@@ -605,9 +606,11 @@ namespace HappyHeart
*
* \warning Only Perform() should call this one!
*/
template<differential is_differential>
void SetIndexInnerLoop(unsigned int index) noexcept;
//! Returns the string that gives current time and current inner loop index.
template<differential is_differential>
std::string GetIterationTag() const noexcept;
public:
......@@ -663,6 +666,11 @@ namespace HappyHeart
//! Non constant accessor to the deltaSolution ( \todo #820 Rename and comment!)
GlobalVector& GetNonCstDeltaSolution() noexcept;
//! Constant accessor to the delta displacement (dispDelta in Freefem).
const GlobalVector& GetDeltaDisplacement() const noexcept;
//! Non constant accessor to the delta displacement (dispDelta in Freefem).
GlobalVector& GetNonCstDeltaDisplacement() noexcept;
private:
......@@ -688,7 +696,8 @@ namespace HappyHeart
//! Compute rhs.
void ComputeRhs(differential is_differential);
template<differential is_differential>
void ComputeRhs();
//! Compute T11.
void ComputeT11(Wrappers::Petsc::DoReuseMatrix do_reuse_matrix,
......@@ -989,6 +998,9 @@ namespace HappyHeart
//! deltaSolution (\todo #820 Rename and comment!) dispDeltaNew in Freefem
GlobalVector::unique_ptr delta_solution_ = nullptr;
//! delta displacement (dispDelta in Freefem).
GlobalVector::unique_ptr delta_displacement_ = nullptr;
///@}
private:
......@@ -1008,8 +1020,8 @@ namespace HappyHeart
//! Internal friction (Dporo in Freefem script).
double internal_friction_ = std::numeric_limits<double>::lowest();
//! Index of the current inner loop (for outputs).
unsigned int index_inner_loop_ = NumericNS::UninitializedIndex<unsigned int>();
//! Index of the current inner loop (for outputs). First is for non differential case, second for differential.
std::array<unsigned int, 2u> index_inner_loop_ {{ NumericNS::UninitializedIndex<unsigned int>(), NumericNS::UninitializedIndex<unsigned int>() }};
//! Parameter that encapsulates velocity part of the current solution.
ParameterAtDof<ParameterNS::Type::vector>::type::const_unique_ptr velocity_solution_param_ = nullptr;
......
......@@ -32,7 +32,7 @@ namespace HappyHeart
double VariationalFormulation<HyperelasticLawT>
::Perform(unsigned int iteration)
{
SetIndexInnerLoop(iteration);
SetIndexInnerLoop<differential::no>(iteration);
decltype(auto) time_manager = parent::GetTimeManager();
decltype(auto) god_of_dof = parent::GetGodOfDof();
decltype(auto) numbering_subset =
......@@ -78,14 +78,14 @@ namespace HappyHeart
parent::GetOutputDirectory(numbering_subset)
+ "/"
+ DifferentialPreffix(differential::no)
+ "full_matrix_" + GetIterationTag() + ".hhdata",
+ "full_matrix_" + GetIterationTag<differential::no>() + ".hhdata",
__FILE__, __LINE__);
system_matrix.View(parent::MpiHappyHeart(),
parent::GetOutputDirectory(numbering_subset)
+ "/"
+ DifferentialPreffix(differential::no)
+ "full_matrix_" + GetIterationTag() + ".m",
+ "full_matrix_" + GetIterationTag<differential::no>() + ".m",
__FILE__, __LINE__,
PETSC_VIEWER_ASCII_MATLAB);
......@@ -93,10 +93,10 @@ namespace HappyHeart
parent::GetOutputDirectory(numbering_subset)
+ "/"
+ DifferentialPreffix(differential::no)
+ "full_matrix_" + GetIterationTag() + ".binary.hhdata",
+ "full_matrix_" + GetIterationTag<differential::no>() + ".binary.hhdata",
__FILE__, __LINE__);
ComputeRhs(differential::no); // must occur after T11, as it needs the t11 matrix.
ComputeRhs<differential::no>(); // must occur after T11, as it needs the t11 matrix.
parent::template SolveLinear<IsFactorized::no>(numbering_subset, numbering_subset);
......@@ -109,7 +109,7 @@ namespace HappyHeart
parent::GetOutputDirectory(numbering_subset)
+ "/"
+ DifferentialPreffix(differential::no)
+ "correction_" + GetIterationTag() + ".m",
+ "correction_" + GetIterationTag<differential::no>() + ".m",
__FILE__, __LINE__,
PETSC_VIEWER_ASCII_MATLAB);
......@@ -127,7 +127,7 @@ namespace HappyHeart
parent::GetOutputDirectory(numbering_subset)
+ "/"
+ DifferentialPreffix(differential::no)
+ "mixt_variables_" + GetIterationTag() + ".m",
+ "mixt_variables_" + GetIterationTag<differential::no>() + ".m",
__FILE__, __LINE__,
PETSC_VIEWER_ASCII_MATLAB);
}
......@@ -148,7 +148,7 @@ namespace HappyHeart
const GlobalMatrix& monolithic_with_dirichlet,
const GlobalMatrix& monolithic_without_dirichlet)
{
SetIndexInnerLoop(iteration);
SetIndexInnerLoop<differential::yes>(iteration);
//decltype(auto) time_manager = parent::GetTimeManager();
decltype(auto) god_of_dof = parent::GetGodOfDof();
decltype(auto) numbering_subset =
......@@ -172,7 +172,7 @@ namespace HappyHeart
const auto filename = parent::GetOutputDirectory(numbering_subset)
+ "/"
+ DifferentialPreffix(differential::yes)
+ "rhs_" + GetIterationTag() + ".m";
+ "rhs_" + GetIterationTag<differential::yes>() + ".m";
rhs.View(parent::MpiHappyHeart(),
filename,
......@@ -202,7 +202,7 @@ namespace HappyHeart
const auto filename = parent::GetOutputDirectory(numbering_subset)
+ "/"
+ DifferentialPreffix(differential::yes)
+ "dirichletFunctionMixtVariables_" + GetIterationTag() + ".m";
+ "dirichletFunctionMixtVariables_" + GetIterationTag<differential::yes>() + ".m";
robin_vector.View(parent::MpiHappyHeart(),
filename,
......@@ -224,7 +224,7 @@ namespace HappyHeart
const auto filename = parent::GetOutputDirectory(numbering_subset)
+ "/"
+ DifferentialPreffix(differential::yes)
+ "dirichletLoadStateCorrection_" + GetIterationTag() + ".m";
+ "dirichletLoadStateCorrection_" + GetIterationTag<differential::yes>() + ".m";
dirichlet_load_state_correction.View(parent::MpiHappyHeart(),
filename,
......@@ -244,7 +244,7 @@ namespace HappyHeart
const auto filename = parent::GetOutputDirectory(numbering_subset)
+ "/"
+ DifferentialPreffix(differential::yes)
+ "imageMixtVariables_just_before_solve_" + GetIterationTag() + ".m";
+ "imageMixtVariables_just_before_solve_" + GetIterationTag<differential::yes>() + ".m";
rhs.View(parent::MpiHappyHeart(),
filename,
......@@ -264,7 +264,7 @@ namespace HappyHeart
const auto filename = parent::GetOutputDirectory(numbering_subset)
+ "/"
+ DifferentialPreffix(differential::yes)
+ "deltaMixtVariables_without_dirichletFunctionMixtVariables_" + GetIterationTag() + ".m";
+ "deltaMixtVariables_without_dirichletFunctionMixtVariables_" + GetIterationTag<differential::yes>() + ".m";
parent::GetSystemSolution(numbering_subset).View(parent::MpiHappyHeart(),
filename,
......@@ -285,7 +285,7 @@ namespace HappyHeart
const auto filename = parent::GetOutputDirectory(numbering_subset)
+ "/"
+ DifferentialPreffix(differential::yes)
+ "deltaMixtVariables_with_dirichletFunctionMixtVariables_" + GetIterationTag() + ".m";
+ "deltaMixtVariables_with_dirichletFunctionMixtVariables_" + GetIterationTag<differential::yes>() + ".m";
delta_mixt_variables.View(parent::MpiHappyHeart(),
filename,
......@@ -319,7 +319,7 @@ namespace HappyHeart
const auto filename = parent::GetOutputDirectory(numbering_subset)
+ "/"
+ DifferentialPreffix(differential::yes)
+ "deltaDarcyVector_" + GetIterationTag() + ".m";
+ "deltaDarcyVector_" + GetIterationTag<differential::yes>() + ".m";
delta_darcy_vector.View(parent::MpiHappyHeart(),
filename,
......@@ -338,7 +338,7 @@ namespace HappyHeart
const auto filename = parent::GetOutputDirectory(numbering_subset)
+ "/"
+ DifferentialPreffix(differential::yes)
+ "deltaSolidVector_" + GetIterationTag() + ".m";
+ "deltaSolidVector_" + GetIterationTag<differential::yes>() + ".m";
delta_solid_vector.View(parent::MpiHappyHeart(),
filename,
......@@ -351,7 +351,7 @@ namespace HappyHeart
const auto filename = parent::GetOutputDirectory(numbering_subset)
+ "/"
+ DifferentialPreffix(differential::yes)
+ "midpoint_position_" + GetIterationTag() + ".m";
+ "midpoint_position_" + GetIterationTag<differential::yes>() + ".m";
midpoint_position.View(parent::MpiHappyHeart(),
filename,
......@@ -400,7 +400,7 @@ namespace HappyHeart
const auto filename = parent::GetOutputDirectory(numbering_subset)
+ "/"
+ DifferentialPreffix(differential::yes)
+ "delta_residual_" + GetIterationTag() + ".m";
+ "delta_residual_" + GetIterationTag<differential::yes>() + ".m";
std::cout << "WRITE FILENAME " << filename << std::endl;
......@@ -432,7 +432,7 @@ namespace HappyHeart
const auto filename = parent::GetOutputDirectory(numbering_subset)
+ "/"
+ DifferentialPreffix(differential::yes)
+ "solid_varf_sol_" + GetIterationTag() + ".m";
+ "solid_varf_sol_" + GetIterationTag<differential::yes>() + ".m";
delta_sol.View(parent::MpiHappyHeart(),
filename,
......@@ -440,12 +440,49 @@ namespace HappyHeart
PETSC_VIEWER_ASCII_MATLAB);
}
decltype(auto) delta_sol = GetDeltaSolution();
auto& delta_displacement = GetNonCstDeltaDisplacement();
{
std::cout << "END INIT DELTA " << GetIterationTag<differential::yes>() << std::endl;
const auto filename = parent::GetOutputDirectory(numbering_subset)
+ "/"
+ DifferentialPreffix(differential::yes)
+ "end_initial_disp_delta_" + GetIterationTag<differential::yes>() + ".m";
delta_displacement.View(parent::MpiHappyHeart(),
filename,
__FILE__, __LINE__,
PETSC_VIEWER_ASCII_MATLAB);
}
// Wrappers::Petsc::AXPY(-1.,
// delta_sol,
// delta_displacement,
// __FILE__, __LINE__);
{
const auto filename = parent::GetOutputDirectory(numbering_subset)
+ "/"
+ DifferentialPreffix(differential::yes)
+ "disp_delta_" + GetIterationTag<differential::yes>() + ".m";
delta_displacement.View(parent::MpiHappyHeart(),
filename,
__FILE__, __LINE__,
PETSC_VIEWER_ASCII_MATLAB);
}
// \todo #820 At the moment, I accept my norm is not the same as Bruno's; see what Petsc gives when
// SNES is applied (system_rhs's is what we got in hyperelastic model).
const double norm = parent::GetSystemRhs(numbering_subset).Norm(NORM_2, __FILE__, __LINE__);
const double norm = delta_displacement.Norm(NORM_2, __FILE__, __LINE__);
return norm;
}
......@@ -551,22 +588,27 @@ namespace HappyHeart
template<class HyperelasticLawT>
template<differential is_differential>
void VariationalFormulation<HyperelasticLawT>::SetIndexInnerLoop(unsigned int index) noexcept
{
index_inner_loop_ = index;
index_inner_loop_[EnumUnderlyingType(is_differential)] = index;
}
template<class HyperelasticLawT>
template<differential is_differential>
std::string VariationalFormulation<HyperelasticLawT>::GetIterationTag() const noexcept
{
return "time_" + std::to_string(parent::GetTimeManager().NtimeModified()) + "_it_"
+ std::to_string(GetIndexInnerLoop());
return "time_" + std::to_string(parent::GetTimeManager().NtimeModified()) + "_"
+ DifferentialPreffix(is_differential)
+ "it_"
+ std::to_string(GetIndexInnerLoop<is_differential>());
}
template<class HyperelasticLawT>
void VariationalFormulation<HyperelasticLawT>::ComputeRhs(differential is_differential)
template<differential is_differential>
void VariationalFormulation<HyperelasticLawT>::ComputeRhs()
{
decltype(auto) god_of_dof = parent::GetGodOfDof();
decltype(auto) bc_manager = DirichletBoundaryConditionManager::GetInstance();
......@@ -584,7 +626,7 @@ namespace HappyHeart
parent::GetOutputDirectory(numbering_subset)
+ "/"
+ DifferentialPreffix(is_differential)
+ "darcy_" + GetIterationTag() + ".m",
+ "darcy_" + GetIterationTag<is_differential>() + ".m",
__FILE__, __LINE__,
PETSC_VIEWER_ASCII_MATLAB);
......@@ -597,7 +639,7 @@ namespace HappyHeart
parent::GetOutputDirectory(numbering_subset)
+ "/"
+ DifferentialPreffix(is_differential)
+ "pressure_contrib_to_rhs_" + GetIterationTag() + ".m",
+ "pressure_contrib_to_rhs_" + GetIterationTag<is_differential>() + ".m",
__FILE__, __LINE__,
PETSC_VIEWER_ASCII_MATLAB);
}
......@@ -651,7 +693,7 @@ namespace HappyHeart
parent::GetOutputDirectory(numbering_subset)
+ "/"
+ DifferentialPreffix(is_differential)
+ "fluid_mass_contrib_to_rhs_" + GetIterationTag() + ".m",
+ "fluid_mass_contrib_to_rhs_" + GetIterationTag<is_differential>() + ".m",
__FILE__, __LINE__,
PETSC_VIEWER_ASCII_MATLAB);
......@@ -674,7 +716,7 @@ namespace HappyHeart
parent::GetOutputDirectory(numbering_subset)
+ "/"
+ DifferentialPreffix(is_differential)
+ "rhs_" + GetIterationTag() + ".m",
+ "rhs_" + GetIterationTag<is_differential>() + ".m",
__FILE__, __LINE__,
PETSC_VIEWER_ASCII_MATLAB);
}
......@@ -740,16 +782,38 @@ namespace HappyHeart
template<class HyperelasticLawT>
void VariationalFormulation<HyperelasticLawT>::SetSolidVelocityForDifferential(Vec petsc_solid_velocity)
{
auto& work_vel = *work_solid_velocity_;
auto& delta_displacement = GetNonCstDeltaDisplacement();
auto& work_vel_fluid = GetNonCstWorkFluidVelocity();
work_vel.SetFromPetscVec(petsc_solid_velocity, __FILE__, __LINE__);
work_vel.Scale(parent::GetTimeManager().GetInverseTimeStep(), __FILE__, __LINE__);
delta_displacement.SetFromPetscVec(petsc_solid_velocity, __FILE__, __LINE__);
{
decltype(auto) god_of_dof = parent::GetGodOfDof();
decltype(auto) numbering_subset =
god_of_dof.GetNumberingSubset(EnumUnderlyingType(NumberingSubsetIndex::fluid_mass_velocity_pressure));
std::cout << "INIT DELTA " << GetIterationTag<differential::yes>() << std::endl;
const auto filename = parent::GetOutputDirectory(numbering_subset)
+ "/"
+ DifferentialPreffix(differential::yes)
+ "initial_disp_delta_" + GetIterationTag<differential::yes>() + ".m";
delta_displacement.View(parent::MpiHappyHeart(),
filename,
__FILE__, __LINE__,
PETSC_VIEWER_ASCII_MATLAB);
}
Wrappers::Petsc::MatMult(variable_holder_parent::GetVariableHolder().GetSolidOnFluidMesh().GetMatrixSolidToFluid(),
work_vel,
delta_displacement,
work_vel_fluid,
__FILE__, __LINE__);
work_vel_fluid.Scale(parent::GetTimeManager().GetInverseTimeStep(), __FILE__, __LINE__);
auto& solid_differential_velocity = GetNonCstSolidDifferentialVelocity();
......
......@@ -480,9 +480,10 @@ namespace HappyHeart
template<class HyperelasticLawT>
template<differential is_differential>
unsigned int VariationalFormulation<HyperelasticLawT>::GetIndexInnerLoop() const noexcept
{
return index_inner_loop_;
return index_inner_loop_[EnumUnderlyingType(is_differential)];
}
......@@ -819,6 +820,24 @@ namespace HappyHeart
return tangent_solid_varf_;
}
template<class HyperelasticLawT>
inline const GlobalVector& VariationalFormulation<HyperelasticLawT>
::GetDeltaDisplacement() const noexcept
{
assert(!(!delta_displacement_));
return *delta_displacement_;
}
template<class HyperelasticLawT>
inline GlobalVector& VariationalFormulation<HyperelasticLawT>
::GetNonCstDeltaDisplacement() noexcept
{
return const_cast<GlobalVector&>(GetDeltaDisplacement());
}
} // namespace InnerLoopNS
......
......@@ -191,7 +191,7 @@ namespace HappyHeart
AllocateGlobalVector(solid_god_of_dof, GetNonCstDeltaResidual());
delta_solution_ = std::make_unique<GlobalVector>(GetDeltaResidual());
delta_displacement_ = std::make_unique<GlobalVector>(GetDeltaResidual());
}
{
......
......@@ -252,11 +252,12 @@ namespace HappyHeart
// We can't take directly the internal pointer of a wrapped vector which pattern would be similar to input_displ_vector
// : residual computation is wrong if we do.
// auto error_code = VecCopy(*HappyHeart vector*, petsc_output_vector);
//
// if (error_code)
// throw Wrappers::Petsc::ExceptionNS::Exception(error_code, "VecCopy",
// __FILE__, __LINE__);
auto error_code = VecCopy(implicit_step_varf.GetInnerLoopVarf().GetDeltaDisplacement().Internal(),
petsc_output_vector);
if (error_code)
throw Wrappers::Petsc::ExceptionNS::Exception(error_code, "VecCopy",
__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