#820 Start introducing a dev method to properly write on disk the global...

#820 Start introducing a dev method to properly write on disk the global linear algebra to be compared with Freefem (it's solely for dev purposes).

Sources/ModelInstances/UnderDevelopment/Poromechanics/ImplicitStepFluid/InnerLoop/VariationalFormulation.hxx

@@ -603,7 +603,7 @@ namespace HappyHeart
                             ret += std::to_string(GetIndexInnerLoop());
                             break;
                         case differential::yes:
-                            ret += std::to_string(variable_holder_parent::GetVariableHolder().gmres_iteration_index);
+                            ret += std::to_string(variable_holder_parent::GetVariableHolder().dh_index);
                     }
                     
                     return ret;

Sources/ModelInstances/UnderDevelopment/Poromechanics/Model.hxx

@@ -34,7 +34,11 @@ namespace HappyHeart
             double absolute_error = 1.e20;
             double relative_error = 1.e20;
             
-            auto iteration_index = 0u;
+            decltype(auto) variable_holder = GetNonCstVariableHolder();
+            auto& iteration_index = variable_holder.implicit_loop_index;
+            iteration_index = 0u;
+            
+            
 
             while (absolute_error > 1.e-10 && relative_error > 1.e-6 && iteration_index < 1u) // \todo #820 Should not remain hardcoded!
                                                                                               // \todo #820 iteration index sup value is obviously temporary for dev purposes...
@@ -52,10 +56,8 @@ namespace HappyHeart
                 
                 UpdateSolidVector();
                 
-                GetVariableHolder().GetFluidMassVector().View(parent::MpiHappyHeart(),
-                                                              parent::GetOutputDirectory() + "/fluidmass_before_solid_solve_time_" + std::to_string(parent::GetTimeManager().NtimeModified()) + ".m",
-                                                              __FILE__, __LINE__,
-                                                              PETSC_VIEWER_ASCII_MATLAB);
+                variable_holder.template DevPrint<DevPhase::none>(variable_holder.GetFluidMassVector(),
+                                                                  "fluidmass_before_solid_solve");
                         
                 SolidVariationalFormulationPolicyT::Solve();
 
@@ -92,7 +94,7 @@ namespace HappyHeart
                 
                 delta_displacement.ZeroEntries(__FILE__, __LINE__);
                 
-                GetNonCstVariableHolder().gmres_iteration_index = static_cast<unsigned int>(-1);
+                GetNonCstVariableHolder().dh_index = static_cast<unsigned int>(-1);
                 dH_solver_->SolveLinear(gmres_shell_matrix, gmres_shell_matrix,
                                         GetVector<SolidIndex::H0>(),
                                         delta_displacement,
@@ -242,7 +244,7 @@ namespace HappyHeart
         void Model<SolidVariationalFormulationPolicyT>
         ::ShellMatrixOperation(Vec petsc_input_vector, Vec& petsc_output_vector)
         {
-            ++GetNonCstVariableHolder().gmres_iteration_index;
+            ++GetNonCstVariableHolder().dh_index;
             
             auto& implicit_step_varf = GetNonCstImplicitStepFluidVarf();
             

Sources/ModelInstances/UnderDevelopment/Poromechanics/ModelInitialize.hxx

@@ -44,7 +44,8 @@ namespace HappyHeart
         void Model<SolidVariationalFormulationPolicyT>::SupplInitialize(const InputParameterList& input_parameter_data)
         {
             decltype(auto) time_manager = parent::GetTimeManager();
-            variable_holder_ = std::make_unique<VariableHolder>(time_manager, input_parameter_data);
+            variable_holder_ = std::make_unique<VariableHolder>(parent::MpiHappyHeart(),
+                                                                time_manager, input_parameter_data);
             decltype(auto) solid_god_of_dof = parent::GetGodOfDof(EnumUnderlyingType(MeshIndex::solid));
             
             decltype(auto) fluid_god_of_dof = parent::GetGodOfDof(EnumUnderlyingType(MeshIndex::fluid));

Sources/ModelInstances/UnderDevelopment/Poromechanics/VariableHolder.cpp

@@ -24,8 +24,11 @@ namespace HappyHeart
     {
         
         
-        VariableHolder::VariableHolder(const TimeManager& time_manager,
+        VariableHolder::VariableHolder(const Wrappers::Mpi& mpi,
+                                       const TimeManager& time_manager,
                                        const InputParameterList& input_parameter_data)
+        : Crtp::HappyHeartMpi<VariableHolder>(mpi),
+        time_manager_(time_manager)
         {
             decltype(auto) god_of_dof_manager = GodOfDofManager::GetInstance();
             decltype(auto) fluid_god_of_dof = god_of_dof_manager.GetGodOfDof(EnumUnderlyingType(MeshIndex::fluid));
@@ -33,6 +36,18 @@ namespace HappyHeart
             decltype(auto) fluid_mesh = fluid_god_of_dof.GetGeometricMeshRegion();
             decltype(auto) solid_mesh = solid_god_of_dof.GetGeometricMeshRegion();
 
+            {
+                namespace ipl = Utilities::InputParameterListNS;
+                output_directory_ = ipl::Extract<InputParameter::Result::OutputDirectory>::Folder(input_parameter_data);
+                
+                output_directory_ += "/Dev/";
+                
+                if (FilesystemNS::Folder::DoExist(output_directory_))
+                    FilesystemNS::Folder::Remove(output_directory_, __FILE__, __LINE__);
+                
+                
+                FilesystemNS::Folder::Create(output_directory_, __FILE__, __LINE__);
+            }
             
             fluid_density_ = InitScalarParameter<InputParameter::Fluid::Density>("Fluid density",
                                                                                  fluid_mesh,

Sources/ModelInstances/UnderDevelopment/Poromechanics/VariableHolder.hpp

@@ -31,6 +31,14 @@ namespace HappyHeart
     {
         
         
+        enum class DevPhase
+        {
+            none,
+            newton_fp,
+            dH
+        };
+        
+        
         /*!
          * \brief The purpose of this struct is to hold quantities that are used at many levels 
          *
@@ -42,6 +50,7 @@ namespace HappyHeart
          * the appropriate VariationalFormulation.
          */
         class VariableHolder
+        : public Crtp::HappyHeartMpi<VariableHolder>
         {
             
         public:
@@ -58,7 +67,8 @@ namespace HappyHeart
             ///@{
             
             //! Constructor.
-            explicit VariableHolder(const TimeManager& time_manager,
+            explicit VariableHolder(const Wrappers::Mpi& mpi,
+                                    const TimeManager& time_manager,
                                     const InputParameterList& input_parameter_data);
             
             //! Destructor.
@@ -85,6 +95,11 @@ namespace HappyHeart
             
             const GlobalVector& GetMidpointPosition() const noexcept;
             
+            
+            template<DevPhase DevPhaseT, class GlobalLinearAlgebraT>
+            void DevPrint(const GlobalLinearAlgebraT& linear_algebra, std::string&& filename) const;
+            
+            
         public:
             
             //! Accessor to mass vector.
@@ -126,8 +141,14 @@ namespace HappyHeart
             const Solid& GetSolid() const noexcept;
             
             
-            //! Index of gmres iteration. \todo #820 For dev purpose!
-            unsigned int gmres_iteration_index = 0u;
+            //! Debug indexes...
+            unsigned int implicit_loop_index = 0u;
+            
+            unsigned int newton_fp_index = 0u;
+            
+            unsigned int dh_index = 0u;
+            
+
 
             
         private:
@@ -155,6 +176,13 @@ namespace HappyHeart
             
             //! Midpoint position. // \todo #820 TMP: midpoint time scheme should not be the sole choice!
             const GlobalVector* midpoint_position_ = nullptr;
+            
+            //! Output directory (for dev).
+            std::string output_directory_;
+            
+            //!
+            const TimeManager& time_manager_;
+        
          
         };
         

Sources/ModelInstances/UnderDevelopment/Poromechanics/VariableHolder.hxx

@@ -98,7 +98,39 @@ namespace HappyHeart
             assert(!(!midpoint_position_) && "Make sure SetMidpointPosition() was correctly called!");
             return *midpoint_position_;
         }
+        
 
+        template<DevPhase DevPhaseT, class GlobalLinearAlgebraT>
+        void VariableHolder::DevPrint(const GlobalLinearAlgebraT& linear_algebra,
+                                   std::string&& filename) const
+        {
+            std::ostringstream oconv;
+            
+            oconv << output_directory_ << filename
+            << "_time_it_" << time_manager_.NtimeModified()
+            << "_implicit_loop_" << implicit_loop_index;
+            
+            switch(DevPhaseT)
+            {
+                case DevPhase::none:
+                    break;
+                case DevPhase::newton_fp:
+                    oconv << "_newton_fp_" << newton_fp_index;
+                    break;
+                case DevPhase::dH:
+                    oconv << "_dH_" << dh_index;
+                    break;
+            }
+            
+            oconv << ".m";
+            
+            std::cout << "FILE = " << oconv.str() << std::endl;
+            
+            linear_algebra.View(MpiHappyHeart(),
+                                oconv.str(),
+                                __FILE__, __LINE__,
+                                PETSC_VIEWER_ASCII_MATLAB);
+        }