#723 FSI/Newton: introduce a stop criterion. I didn't use exactly the same as...

#723 FSI/Newton: introduce a stop criterion. I didn't use exactly the same as freefem, but it should do (especially considering I'll use SNES later on to handle the whole coupling and let it care about that).

Data/Lua/demo_input_fsi_ei.lua

@@ -21,7 +21,7 @@ transient = {
 	-- Maximum time, if set to zero run a case.
 	-- Expected format: VALUE
 	-- Constraint: v >= 0.
-	timeMax = 0.5e-4
+	timeMax = 3.1e-4
 
 } -- transient
 

Sources/ModelInstances/FSI_EI/Newton/Model.hpp

@@ -167,7 +167,7 @@ namespace HappyHeart
              * performed within the method because differential and non differential cases use different values
              * for this velocity.
              */
-            void NewUpdateFluidDirichletCondition();
+            void UpdateFluidDirichletCondition();
             
             //! Compute solid residual (in input of solid variational formulation).
             void ComputeSolidResidual();

Sources/ModelInstances/FSI_EI/Newton/Model.hxx

@@ -308,6 +308,11 @@ namespace HappyHeart
                                          __FILE__, __LINE__);
             }
             
+            {
+                GetNonCstVector<Solid::dispSOldr_on_interface>().Copy(GetVector<Solid::dispSr_on_interface>(),
+                                                                      __FILE__, __LINE__);
+            }
+            
          
         }
         

Sources/ModelInstances/FSI_EI/Newton/ModelForward.hxx

@@ -38,8 +38,13 @@ namespace HappyHeart
             ExplicitFluidStep(output_dir);
             
             const auto& gmres_shell_matrix = GetGmresShellMatrix();
+            
+            double absolute_error = 1.e20;
+            double relative_error = 1.e20;
+            
+            auto iteration_index = 0u;
 
-            for (unsigned int i = 0u; i < 3u; ++i) // \todo #723 Implement stop conditions! Or even better define a SNES for that...
+            while (absolute_error > 1.e-10 && relative_error > 1.e-6 && iteration_index < 1000u) // \todo #723 Should not remain hardcoded!
             {
                
                 GetNonCstVector<Solid::dispSrk0>().Copy(GetVector<Solid::dispSr>(), __FILE__, __LINE__);
@@ -64,7 +69,7 @@ namespace HappyHeart
                     velFromSr_on_interface.Scale(inv_time_step, __FILE__, __LINE__);
                 }
                 
-                NewUpdateFluidDirichletCondition();
+                UpdateFluidDirichletCondition();
                 
                 ImplicitFluidStep(differential::no);
                 
@@ -135,15 +140,28 @@ namespace HappyHeart
                                           __FILE__, __LINE__);
                 }
                 
+                const double norm = Wrappers::Petsc::NormL2(mpi, dispSdelta, __FILE__, __LINE__);
+                
+                if (iteration_index == 0u)
+                    absolute_error = norm;
+                else
+                {
+                    assert(absolute_error > 1.e-10); // \todo #723 Should be here exactly same condition as in the while.
+                    
+                    relative_error = norm / absolute_error;
+                }
                 
+                std::cout << "End of coupling loop "<< iteration_index << " - abs = "<< absolute_error << " and "
+                "rel = " << relative_error << std::endl;
                 
+                ++iteration_index;
                 
             } // Coupling Newton loop
         }
         
         
         template<class SolidVariationalFormulationPolicyT>
-        void Model<SolidVariationalFormulationPolicyT>::NewUpdateFluidDirichletCondition()
+        void Model<SolidVariationalFormulationPolicyT>::UpdateFluidDirichletCondition()
         {
             const auto& velFromSr_on_interface = GetVector<Fluid::velFromSr_on_interface>();
             auto& p2_solid_velocity_on_interface = GetNonCstVector<Fluid::p2_vel_on_interface>();
@@ -189,7 +207,7 @@ namespace HappyHeart
                                                             __FILE__, __LINE__);
                     
                     velFromSr_on_interface.Copy(dispSr_on_interface, __FILE__, __LINE__);
-                    Wrappers::Petsc::AXPY(-1.,  // \todo #723  NOT SURE AT ALL IF THIS PART SHOULD BE INTRODUCED here; if so NewUpdateFluidDirichletCondition() can be improved (take input displacement in parameters).
+                    Wrappers::Petsc::AXPY(-1.,  // \todo #723  NOT SURE AT ALL IF THIS PART SHOULD BE INTRODUCED here; if so UpdateFluidDirichletCondition() can be improved (take input displacement in parameters).
                                           GetVector<Solid::dispSOldr_on_interface>(),
                                           velFromSr_on_interface,
                                           __FILE__, __LINE__);
@@ -199,7 +217,7 @@ namespace HappyHeart
                 }
 
                 
-                NewUpdateFluidDirichletCondition();
+                UpdateFluidDirichletCondition();
             }
             
             ImplicitFluidStep(differential::yes);

Sources/ThirdParty/Wrappers/Petsc/Vector/Vector.cpp

@@ -503,6 +503,25 @@ namespace HappyHeart
             }
             
             
+            double NormL2(const Mpi& mpi,
+                          const Vector& v,
+                          const char* invoking_file, int invoking_line)
+            {
+                AccessVectorContent<Utilities::Access::read_only> content(v, invoking_file, invoking_line);
+                
+                const auto size = content.GetSize(invoking_file, invoking_line);
+                
+                double ret = 0.;
+                
+                for (auto i = 0u; i < size; ++i)
+                    ret += Utilities::Square(content.GetValue(i));
+                
+                // Query all processors
+                return std::sqrt(mpi.AllReduce(ret, MpiNS::Op::Sum));
+            }
+
+            
+            
             void Vector::SetFromPetscVec(const Vec& petsc_vector,
                                          const char* invoking_file, int invoking_line)
             {

Sources/ThirdParty/Wrappers/Petsc/Vector/Vector.hpp

@@ -598,6 +598,12 @@ namespace HappyHeart
                                  const Vector& v2,
                                  const char* invoking_file, int invoking_line);
             
+            /*!
+             * \brief Compute Norm L2.
+             */
+            double NormL2(const Mpi& mpi,
+                          const Vector& vector,
+                          const char* invoking_file, int invoking_line);
 
    
             template<Utilities::Access AccessT>