#820 In implicit step inner loop, introduce a method to handle the...

#820 In implicit step inner loop, introduce a method to handle the differential case, for which several operations are very different of the non differential case (some operators are for instance only a limited subpart of the non differential case).

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;

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>

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!
                 

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.