diff --git a/Sources/ModelInstances/UnderDevelopment/Poromechanics/ImplicitStepFluid/InnerLoop/VariationalFormulation.hpp b/Sources/ModelInstances/UnderDevelopment/Poromechanics/ImplicitStepFluid/InnerLoop/VariationalFormulation.hpp
index 0902d60377b42d022274d61083dd4072f177ec92..3c35804f1d25d626d78e295357216f20dc277b82 100644
--- a/Sources/ModelInstances/UnderDevelopment/Poromechanics/ImplicitStepFluid/InnerLoop/VariationalFormulation.hpp
+++ b/Sources/ModelInstances/UnderDevelopment/Poromechanics/ImplicitStepFluid/InnerLoop/VariationalFormulation.hpp
@@ -174,7 +174,15 @@ namespace HappyHeart
*
* \return Norm obtained.
*/
- double Perform(differential is_differential, unsigned int iteration);
+ double Perform(unsigned int iteration);
+
+
+ /*!
+ *
+ * \todo #9 - #820 To complete!
+ */
+ double PerformDifferential(unsigned int iteration);
+
//! Compute the matrix with an updated T21 block (to use after the Newton converged).
const GlobalMatrix& ComputeMatrixWithUpdatedT21() noexcept;
diff --git a/Sources/ModelInstances/UnderDevelopment/Poromechanics/ImplicitStepFluid/InnerLoop/VariationalFormulation.hxx b/Sources/ModelInstances/UnderDevelopment/Poromechanics/ImplicitStepFluid/InnerLoop/VariationalFormulation.hxx
index d35e5cee71cbf5da569177e361df457eaa91795d..9c233da2ecb02a579ab1a54e0f57b2d60c66436e 100644
--- a/Sources/ModelInstances/UnderDevelopment/Poromechanics/ImplicitStepFluid/InnerLoop/VariationalFormulation.hxx
+++ b/Sources/ModelInstances/UnderDevelopment/Poromechanics/ImplicitStepFluid/InnerLoop/VariationalFormulation.hxx
@@ -30,8 +30,7 @@ namespace HappyHeart
template<class HyperelasticLawT>
double VariationalFormulation<HyperelasticLawT>
- ::Perform(differential is_differential,
- unsigned int iteration)
+ ::Perform(unsigned int iteration)
{
SetIndexInnerLoop(iteration);
decltype(auto) time_manager = parent::GetTimeManager();
@@ -68,49 +67,36 @@ namespace HappyHeart
corrected_values.UpdateGhosts(__FILE__, __LINE__);
}
- switch(is_differential)
- {
- case differential::no:
- {
- auto& system_matrix = parent::GetNonCstSystemMatrix(numbering_subset, numbering_subset);
- system_matrix.ZeroEntries(__FILE__, __LINE__);
-
- ComputeMatrix(do_reuse_matrix, BlocksToRecompute::all);
-
- system_matrix.Assembly(__FILE__, __LINE__);
-
- system_matrix.View(parent::MpiHappyHeart(),
- parent::GetOutputDirectory(numbering_subset)
- + "/"
- + DifferentialPreffix(is_differential)
- + "full_matrix_" + GetIterationTag() + ".hhdata",
- __FILE__, __LINE__);
-
- system_matrix.View(parent::MpiHappyHeart(),
- parent::GetOutputDirectory(numbering_subset)
- + "/"
- + DifferentialPreffix(is_differential)
- + "full_matrix_" + GetIterationTag() + ".m",
- __FILE__, __LINE__,
- PETSC_VIEWER_ASCII_MATLAB);
-
- system_matrix.ViewBinary(parent::MpiHappyHeart(),
- parent::GetOutputDirectory(numbering_subset)
- + "/"
- + DifferentialPreffix(is_differential)
- + "full_matrix_" + GetIterationTag() + ".binary.hhdata",
- __FILE__, __LINE__);
-
-
- break;
- }
- case differential::yes:
- // Just reuse tangent matrix computed at the end of last non differential case.
- break;
- } // switch
+ auto& system_matrix = parent::GetNonCstSystemMatrix(numbering_subset, numbering_subset);
+ system_matrix.ZeroEntries(__FILE__, __LINE__);
-
- ComputeRhs(is_differential); // must occur after T11, as it needs the t11 matrix.
+ ComputeMatrix(do_reuse_matrix, BlocksToRecompute::all);
+
+ system_matrix.Assembly(__FILE__, __LINE__);
+
+ system_matrix.View(parent::MpiHappyHeart(),
+ parent::GetOutputDirectory(numbering_subset)
+ + "/"
+ + DifferentialPreffix(differential::no)
+ + "full_matrix_" + GetIterationTag() + ".hhdata",
+ __FILE__, __LINE__);
+
+ system_matrix.View(parent::MpiHappyHeart(),
+ parent::GetOutputDirectory(numbering_subset)
+ + "/"
+ + DifferentialPreffix(differential::no)
+ + "full_matrix_" + GetIterationTag() + ".m",
+ __FILE__, __LINE__,
+ PETSC_VIEWER_ASCII_MATLAB);
+
+ system_matrix.ViewBinary(parent::MpiHappyHeart(),
+ parent::GetOutputDirectory(numbering_subset)
+ + "/"
+ + DifferentialPreffix(differential::no)
+ + "full_matrix_" + GetIterationTag() + ".binary.hhdata",
+ __FILE__, __LINE__);
+
+ ComputeRhs(differential::no); // must occur after T11, as it needs the t11 matrix.
parent::template SolveLinear<IsFactorized::no>(numbering_subset, numbering_subset);
@@ -122,7 +108,7 @@ namespace HappyHeart
solution.View(parent::MpiHappyHeart(),
parent::GetOutputDirectory(numbering_subset)
+ "/"
- + DifferentialPreffix(is_differential)
+ + DifferentialPreffix(differential::no)
+ "correction_" + GetIterationTag() + ".m",
__FILE__, __LINE__,
PETSC_VIEWER_ASCII_MATLAB);
@@ -140,7 +126,7 @@ namespace HappyHeart
corrected_values.View(parent::MpiHappyHeart(),
parent::GetOutputDirectory(numbering_subset)
+ "/"
- + DifferentialPreffix(is_differential)
+ + DifferentialPreffix(differential::no)
+ "mixt_variables_" + GetIterationTag() + ".m",
__FILE__, __LINE__,
PETSC_VIEWER_ASCII_MATLAB);
@@ -156,6 +142,22 @@ namespace HappyHeart
}
+ template<class HyperelasticLawT>
+ double VariationalFormulation<HyperelasticLawT>
+ ::PerformDifferential(unsigned int iteration)
+ {
+ SetIndexInnerLoop(iteration);
+ decltype(auto) time_manager = parent::GetTimeManager();
+ decltype(auto) god_of_dof = parent::GetGodOfDof();
+ decltype(auto) numbering_subset =
+ god_of_dof.GetNumberingSubset(EnumUnderlyingType(NumberingSubsetIndex::fluid_mass_velocity_pressure));
+
+ // \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__);
+ return norm;
+ }
+
template<class HyperelasticLawT>
void VariationalFormulation<HyperelasticLawT>
diff --git a/Sources/ModelInstances/UnderDevelopment/Poromechanics/ImplicitStepFluid/VariationalFormulation.hxx b/Sources/ModelInstances/UnderDevelopment/Poromechanics/ImplicitStepFluid/VariationalFormulation.hxx
index 7df9caf3b335ec7fb171ec62f0476a9bf7951de8..2b010e3310b52c7c3c48ef8bdbd62401fafcfb23 100644
--- a/Sources/ModelInstances/UnderDevelopment/Poromechanics/ImplicitStepFluid/VariationalFormulation.hxx
+++ b/Sources/ModelInstances/UnderDevelopment/Poromechanics/ImplicitStepFluid/VariationalFormulation.hxx
@@ -72,7 +72,20 @@ namespace HappyHeart
// Newton loop (\todo #820 To replace with Petsc Newton...)
do
{
- current_norm = inner_varf.Perform(is_differential, counter++);
+ switch(is_differential)
+ {
+ case differential::no:
+ {
+ current_norm = inner_varf.Perform(counter++);
+ break;
+ }
+ case differential::yes:
+ {
+ current_norm = inner_varf.PerformDifferential(counter++);
+ break;
+ }
+ }
+
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!
diff --git a/Sources/ModelInstances/UnderDevelopment/Poromechanics/Model.hxx b/Sources/ModelInstances/UnderDevelopment/Poromechanics/Model.hxx
index 2d2c75c2cc777dc9ce40d87d9e8ba940e3d8c2bb..1ebe354177d78d26bdfe5a364820a8a0adc22f32 100644
--- a/Sources/ModelInstances/UnderDevelopment/Poromechanics/Model.hxx
+++ b/Sources/ModelInstances/UnderDevelopment/Poromechanics/Model.hxx
@@ -247,6 +247,10 @@ namespace HappyHeart
input_vector.Scale(parent::GetTimeManager().GetInverseTimeStep(), __FILE__, __LINE__);
+ auto& implicit_step_varf = GetNonCstImplicitStepFluidVarf();
+
+ implicit_step_varf.Perform(differential::yes);
+
// 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.