#902 Add safety whenever MatAXPY is called to check whether same non zero...

#902 Add safety whenever MatAXPY is called to check whether same non zero pattern might be used. The default value has been disabled: this parameter is way too important to be hidden.

Data/Lua/demo_input_fsi_ei_two_meshes.lua

@@ -1125,6 +1125,11 @@ Petsc2 = {
     -- Constraint: v > 0.
     relativeTolerance = 1.e-4,
     
+    -- Step size tolerance
+    -- Expected format: {VALUE1, VALUE2, ...}
+    -- Constraint: v > 0.
+    stepSizeTolerance = 1e-8,
+    
     -- List of solver: { chebychev cg gmres preonly bicg python };
     -- Expected format: {"VALUE1", "VALUE2", ...}
     -- Constraint: ops_in(v, {'Mumps', 'Umfpack', 'Gmres'})

Sources/Core/LinearAlgebra/GlobalMatrix.cpp

@@ -105,6 +105,16 @@ namespace HappyHeart
             ExceptionNS::PrintAndAbort(mpi, e.what());
         }
     }
+    
+    
+    void AssertSameNumberingSubset(const GlobalMatrix& matrix1,
+                                   const GlobalMatrix& matrix2)
+    {
+        assert(matrix1.GetRowNumberingSubset() == matrix2.GetRowNumberingSubset());
+        assert(matrix1.GetColNumberingSubset() == matrix2.GetColNumberingSubset());
+    }
+
+    
     #endif // NDEBUG
 
     

Sources/Core/LinearAlgebra/GlobalMatrix.hpp

@@ -132,6 +132,18 @@ namespace HappyHeart
                                     const GlobalMatrix& matrix,
                                     const Wrappers::Petsc::MatrixPattern& pattern,
                                     const char* invoking_file, int invoking_line);
+    
+    /*!
+     * \brief Assert two matrices share the same \a NumberingSubset.
+     *
+     * \param[in] matrix1 First matrix.
+     * \param[in] matrix2 Second matrix.
+     */
+    void AssertSameNumberingSubset(const GlobalMatrix& matrix1,
+                                   const GlobalMatrix& matrix2);
+    
+    
+    
     # endif // NDEBUG
     
 

Sources/Core/NumberingSubset/NumberingSubset.hpp

@@ -148,7 +148,7 @@ namespace HappyHeart
     //! Comparison use the underlying id.
     bool operator<(const NumberingSubset& lhs, const NumberingSubset& rhs);
     
-    
+   
     ///@} // \addtogroup
 
         

Sources/FiniteElement/FiniteElementSpace/GodOfDof.cpp

@@ -1120,12 +1120,12 @@ namespace HappyHeart
         decltype(auto) pattern = god_of_dof.GetMatrixPattern(matrix.GetRowNumberingSubset(),
                                                              matrix.GetColNumberingSubset());
         
-# ifndef NDEBUG
+        # ifndef NDEBUG
         AssertMatrixRespectPattern(god_of_dof.MpiHappyHeart(),
                                    matrix,
                                    pattern,
                                    invoking_file, invoking_line);
-# endif // NDEBUG
+        # endif // NDEBUG
         
         static_cast<void>(god_of_dof);
         static_cast<void>(matrix);

Sources/FormulationSolver/DofSourcePolicy/DofSource.cpp

@@ -38,7 +38,6 @@ namespace HappyHeart
                 Wrappers::Petsc::AXPY(factor_, dof_source_, rhs, __FILE__, __LINE__);
             }
             
-          
             
         } //namespace DofSourcePolicyNS
         

Sources/Geometry/Mesh/GeometricMeshRegion.cpp

@@ -142,7 +142,9 @@ namespace HappyHeart
     
     
     
-    unsigned int GeometricMeshRegion::Read(const std::string& mesh_file, Internal::MeshNS::FormatNS::Type format, const double space_unit)
+    unsigned int GeometricMeshRegion::Read(const std::string& mesh_file,
+                                           Internal::MeshNS::FormatNS::Type format,
+                                           const double space_unit)
     {
         GeometricElt::vector_shared_ptr unsort_element_list;
         

Sources/ModelInstances/AcousticWave/VariationalFormulation.cpp

@@ -215,7 +215,12 @@ namespace HappyHeart
             system_matrix.Copy(mass_matrix, __FILE__, __LINE__);
             
             const double theta = GetTheta();
-            Wrappers::Petsc::AXPY(theta, stiffness_matrix, system_matrix, __FILE__, __LINE__);
+            
+            #ifndef NDEBUG
+            AssertSameNumberingSubset(stiffness_matrix, system_matrix);
+            #endif // NDEBUG
+            
+            Wrappers::Petsc::AXPY(SAME_NONZERO_PATTERN, theta, stiffness_matrix, system_matrix, __FILE__, __LINE__);
         }
         
         

Sources/ModelInstances/Bidomain/BidomainVariationalFormulation.hxx

@@ -337,13 +337,18 @@ namespace HappyHeart
             const auto& matrix_bidomain = GetMatrixBidomain();
             const auto& matrix_capacity = GetMatrixCapacityPerTimeStep();
 
-            Wrappers::Petsc::AXPY(1., matrix_bidomain, system_matrix, __FILE__, __LINE__);
-            Wrappers::Petsc::AXPY(1., matrix_capacity, system_matrix, __FILE__, __LINE__);
-            
             // Test on the inversibility of the system matrix by adding a mass matrix on Ue.
             auto& mass_matrix = GetMatrixExtracellularPotentialCapacity();
             
-            Wrappers::Petsc::AXPY(1.e-06, mass_matrix, system_matrix, __FILE__, __LINE__);
+            #ifndef NDEBUG
+            AssertSameNumberingSubset(system_matrix, matrix_bidomain);
+            AssertSameNumberingSubset(system_matrix, matrix_capacity);
+            AssertSameNumberingSubset(system_matrix, mass_matrix);
+            #endif // NDEBUG
+            
+            Wrappers::Petsc::AXPY(SAME_NONZERO_PATTERN, 1., matrix_bidomain, system_matrix, __FILE__, __LINE__);
+            Wrappers::Petsc::AXPY(SAME_NONZERO_PATTERN, 1., matrix_capacity, system_matrix, __FILE__, __LINE__);            
+            Wrappers::Petsc::AXPY(SAME_NONZERO_PATTERN, 1.e-06, mass_matrix, system_matrix, __FILE__, __LINE__);
         }
         
         

Sources/ModelInstances/BidomainVentricles/BidomainVentriclesVariationalFormulation.hxx

@@ -352,9 +352,14 @@ namespace HappyHeart
             
             const auto& matrix_bidomain = GetMatrixBidomain();
             const auto& matrix_capacity = GetMatrixCapacityPerTimeStep();
-
-            Wrappers::Petsc::AXPY(1., matrix_bidomain, system_matrix, __FILE__, __LINE__);
-            Wrappers::Petsc::AXPY(1., matrix_capacity, system_matrix, __FILE__, __LINE__);
+            
+            #ifndef NDEBUG
+            AssertSameNumberingSubset(system_matrix, matrix_bidomain);
+            AssertSameNumberingSubset(system_matrix, matrix_capacity);
+            #endif // NDEBUG
+            
+            Wrappers::Petsc::AXPY(SAME_NONZERO_PATTERN, 1., matrix_bidomain, system_matrix, __FILE__, __LINE__);
+            Wrappers::Petsc::AXPY(SAME_NONZERO_PATTERN, 1., matrix_capacity, system_matrix, __FILE__, __LINE__);
         }
         
         
@@ -390,7 +395,7 @@ namespace HappyHeart
             
             rhs.ZeroEntries(__FILE__, __LINE__);
             
-            Wrappers::Petsc::AXPY(1., GetVectorTranscellularIonSource(), rhs, __FILE__, __LINE__);
+            Wrappers::Petsc::AXPY(SAME_NONZERO_PATTERN, 1., GetVectorTranscellularIonSource(), rhs, __FILE__, __LINE__);
         }
         
         
@@ -407,7 +412,7 @@ namespace HappyHeart
                 
                 if (time > time_init && time < time_final)
                 {
-                    Wrappers::Petsc::AXPY(1., GetVectorTranscellularCurrentAppliedOnLeftVentricle(), rhs, __FILE__, __LINE__);
+                    Wrappers::Petsc::AXPY(SAME_NONZERO_PATTERN, 1., GetVectorTranscellularCurrentAppliedOnLeftVentricle(), rhs, __FILE__, __LINE__);
                 }
             }
             
@@ -421,7 +426,7 @@ namespace HappyHeart
                 
                 if (time > time_init && time < time_final)
                 {
-                    Wrappers::Petsc::AXPY(1., GetVectorTranscellularCurrentAppliedOnRightVentricle(), rhs, __FILE__, __LINE__);
+                    Wrappers::Petsc::AXPY(SAME_NONZERO_PATTERN, 1., GetVectorTranscellularCurrentAppliedOnRightVentricle(), rhs, __FILE__, __LINE__);
                 }
             }
         }
@@ -461,7 +466,7 @@ namespace HappyHeart
             
             auto& vector_of_one = GetNonCstVectorOfOneOnExtracellularPotential();
 
-            Wrappers::Petsc::AXPY(-mean, vector_of_one, system_solution, __FILE__, __LINE__);
+            Wrappers::Petsc::AXPY(SAME_NONZERO_PATTERN, -mean, vector_of_one, system_solution, __FILE__, __LINE__);
         }
         
         

Sources/ModelInstances/CardiacMechanicsPrestress/VariationalFormulation.cpp

@@ -692,12 +692,17 @@ namespace HappyHeart
             {
                 coefficient = 0.5;
                 
-                Wrappers::Petsc::AXPY(1.,
+                #ifndef NDEBUG
+                AssertSameNumberingSubset(GetMatrixTangentStiffness(), system_matrix);
+                AssertSameNumberingSubset(GetMatrixMassPerSquareTimeStep(), system_matrix);
+                #endif // NDEBUG
+                
+                Wrappers::Petsc::AXPY(SAME_NONZERO_PATTERN, 1.,
                                       GetMatrixMassPerSquareTimeStep(),
                                       system_matrix, __FILE__, __LINE__);
             }
             
-            Wrappers::Petsc::AXPY(coefficient,
+            Wrappers::Petsc::AXPY(SAME_NONZERO_PATTERN, coefficient,
                                   GetMatrixTangentStiffness(),
                                   system_matrix, __FILE__, __LINE__);
             
@@ -705,7 +710,11 @@ namespace HappyHeart
             
             if (run_case == StaticOrDynamic::static_)
             {
-                Wrappers::Petsc::AXPY(coefficient,
+                #ifndef NDEBUG
+                AssertSameNumberingSubset(GetMatrixTangentFollowingPressure(), system_matrix);
+                #endif // NDEBUG
+                
+                Wrappers::Petsc::AXPY(SAME_NONZERO_PATTERN, coefficient,
                                       GetMatrixTangentFollowingPressure(),
                                       system_matrix, __FILE__, __LINE__);                
             }

Sources/ModelInstances/Elasticity/VariationalFormulationElasticity.cpp

@@ -190,7 +190,11 @@ namespace HappyHeart
                 const auto coefficient =
                     2. * GetVolumicMass().GetConstantValue() / NumericNS::Square(parent::GetTimeManager().GetTimeStep());
             
-                Wrappers::Petsc::AXPY(coefficient,
+                #ifndef NDEBUG
+                AssertSameNumberingSubset(mass, system_matrix);
+                #endif // NDEBUG
+                
+                Wrappers::Petsc::AXPY(SAME_NONZERO_PATTERN, coefficient,
                                       mass,
                                       system_matrix,
                                       __FILE__, __LINE__);
@@ -206,7 +210,11 @@ namespace HappyHeart
                 
                 displacement_previous_time_iteration_matrix.Scale(coefficient, __FILE__, __LINE__);
                 
-                Wrappers::Petsc::AXPY(-.5,
+                #ifndef NDEBUG
+                AssertSameNumberingSubset(stiffness, displacement_previous_time_iteration_matrix);
+                #endif // NDEBUG
+                
+                Wrappers::Petsc::AXPY(SAME_NONZERO_PATTERN, -.5,
                                       stiffness,
                                       displacement_previous_time_iteration_matrix,
                                       __FILE__, __LINE__);

Sources/ModelInstances/FSI_EI/ExplicitStepVariationalFormulation.cpp

@@ -215,7 +215,11 @@ namespace HappyHeart
                                                    GetFluidDensity().GetConstantValue() / time_step);
                 GetMass().Assemble(std::make_tuple(std::ref(matrix)));
                 
-                Wrappers::Petsc::AXPY(1., mass_fluid_velocity, residual_matrix, __FILE__, __LINE__);
+                #ifndef NDEBUG
+                AssertSameNumberingSubset(residual_matrix, mass_fluid_velocity);
+                #endif // NDEBUG
+
+                Wrappers::Petsc::AXPY(SAME_NONZERO_PATTERN, 1., mass_fluid_velocity, residual_matrix, __FILE__, __LINE__);
             }
             
             

Sources/ModelInstances/FSI_EI/SolidVariationalFormulationPolicy/Elasticity/VariationalFormulation.cpp

@@ -121,7 +121,11 @@ namespace HappyHeart
                 const auto coefficient =
                 2. * GetVolumicMass().GetConstantValue() / NumericNS::Square(parent::GetTimeManager().GetTimeStep());
                 
-                Wrappers::Petsc::AXPY(coefficient,
+                #ifndef NDEBUG
+                AssertSameNumberingSubset(mass, system_matrix);
+                #endif // NDEBUG
+                
+                Wrappers::Petsc::AXPY(SAME_NONZERO_PATTERN, coefficient,
                                       mass,
                                       system_matrix,
                                       __FILE__, __LINE__);
@@ -137,7 +141,11 @@ namespace HappyHeart
                 
                 displacement_previous_time_iteration_matrix.Scale(coefficient, __FILE__, __LINE__);
                 
-                Wrappers::Petsc::AXPY(-.5,
+                #ifndef NDEBUG
+                AssertSameNumberingSubset(stiffness, displacement_previous_time_iteration_matrix);
+                #endif // NDEBUG
+                
+                Wrappers::Petsc::AXPY(SAME_NONZERO_PATTERN, -.5,
                                       stiffness,
                                       displacement_previous_time_iteration_matrix,
                                       __FILE__, __LINE__);

Sources/ModelInstances/Heat/HeatVariationalFormulation.cpp

@@ -235,7 +235,10 @@ namespace HappyHeart
                 GetConductivityOperator().Assemble(std::make_tuple(std::ref(matrix)));
             }
             
-            Wrappers::Petsc::AXPY(1., GetMatrixRobin(), system_matrix, __FILE__, __LINE__);
+            #ifndef NDEBUG
+            AssertSameNumberingSubset(GetMatrixRobin(), system_matrix);
+            #endif // NDEBUG
+            Wrappers::Petsc::AXPY(SAME_NONZERO_PATTERN, 1., GetMatrixRobin(), system_matrix, __FILE__, __LINE__);
         }
         
         
@@ -319,8 +322,11 @@ namespace HappyHeart
             const GlobalMatrix& capacity_matrix = GetMatrixCapacityPerTimeStep();
            
             {
-                // System matrix.
-                Wrappers::Petsc::AXPY(1.,
+                #ifndef NDEBUG
+                AssertSameNumberingSubset(capacity_matrix, system_matrix);
+                #endif // NDEBUG
+                
+                Wrappers::Petsc::AXPY(SAME_NONZERO_PATTERN, 1.,
                                       capacity_matrix,
                                       system_matrix,
                                       __FILE__, __LINE__);

Sources/ModelInstances/Hyperelasticity/TimeSchemes/HalfSum/HalfSum.cpp

@@ -30,7 +30,6 @@ namespace HappyHeart
             {
                 static_cast<void>(is_dynamic_iteration);
                 static_cast<void>(vm);
-                
             }
 
             
@@ -62,11 +61,18 @@ namespace HappyHeart
             void ComputeDynamicTimeSchemeTangentContribution(const VectorsAndMatrices<HyperelasticityNS::TimeScheme::half_sum>& vm,
                                                              GlobalMatrix& system_matrix)
             {
-                Wrappers::Petsc::AXPY(0.5,
+                #ifndef NDEBUG
+                AssertSameNumberingSubset(system_matrix, vm.GetMatrixNewStiffness());
+                AssertSameNumberingSubset(system_matrix, vm.GetMatrixStiffnessPreviousTimeIteration());
+                #endif // NDEBUG
+                
+                Wrappers::Petsc::AXPY(SAME_NONZERO_PATTERN,
+                                      0.5,
                                       vm.GetMatrixNewStiffness(),
                                       system_matrix, __FILE__, __LINE__);
                 
-                Wrappers::Petsc::AXPY(0.5,
+                Wrappers::Petsc::AXPY(SAME_NONZERO_PATTERN,
+                                      0.5,
                                       vm.GetMatrixStiffnessPreviousTimeIteration(),
                                       system_matrix, __FILE__, __LINE__);
             }

Sources/ModelInstances/Hyperelasticity/TimeSchemes/Midpoint/Midpoint.cpp

@@ -65,9 +65,13 @@ namespace HappyHeart
 
             
             void ComputeDynamicTimeSchemeTangentContribution(const VectorsAndMatrices<HyperelasticityNS::TimeScheme::midpoint>& vm,
-                                                      GlobalMatrix& system_matrix)
+                                                             GlobalMatrix& system_matrix)
             {
-                Wrappers::Petsc::AXPY(0.5,
+                #ifndef NDEBUG
+                AssertSameNumberingSubset(vm.GetMatrixNewStiffness(), system_matrix);
+                #endif // NDEBUG
+                
+                Wrappers::Petsc::AXPY(SAME_NONZERO_PATTERN, 0.5,
                                       vm.GetMatrixNewStiffness(),
                                       system_matrix, __FILE__, __LINE__);
             }

Sources/ModelInstances/Legacy/FSI_EI_One_Mesh/ExplicitStepVariationalFormulation.cpp

@@ -181,7 +181,10 @@ namespace HappyHeart
                                                    GetFluidDensity().GetConstantValue() / time_step);
                 GetMass().Assemble(std::make_tuple(std::ref(matrix)));
                 
-                Wrappers::Petsc::AXPY(1., mass_fluid_velocity, residual_matrix, __FILE__, __LINE__);
+                #ifndef NDEBUG
+                AssertSameNumberingSubset(mass_fluid_velocity, residual_matrix);
+                #endif // NDEBUG
+                Wrappers::Petsc::AXPY(SAME_NONZERO_PATTERN, 1., mass_fluid_velocity, residual_matrix, __FILE__, __LINE__);
             }
             
             
@@ -226,7 +229,7 @@ namespace HappyHeart
                                      residual,
                                      __FILE__, __LINE__);
             
-            Wrappers::Petsc::AXPY(-1., GetResidualRhs(), residual,
+            Wrappers::Petsc::AXPY(SAME_NONZERO_PATTERN, -1., GetResidualRhs(), residual,
                                   __FILE__, __LINE__);
             
             return residual;

Sources/ModelInstances/Legacy/FSI_EI_One_Mesh/SolidVariationalFormulationPolicy/Elasticity/VariationalFormulation.cpp

@@ -120,7 +120,12 @@ namespace HappyHeart
                 const auto coefficient =
                 2. * GetVolumicMass().GetConstantValue() / NumericNS::Square(parent::GetTimeManager().GetTimeStep());
                 
-                Wrappers::Petsc::AXPY(coefficient,
+                #ifndef NDEBUG
+                AssertSameNumberingSubset(system_matrix, mass);
+                #endif // NDEBUG
+                
+                Wrappers::Petsc::AXPY(SAME_NONZERO_PATTERN,
+                                      coefficient,
                                       mass,
                                       system_matrix,
                                       __FILE__, __LINE__);
@@ -136,7 +141,12 @@ namespace HappyHeart
                 
                 displacement_previous_time_iteration_matrix.Scale(coefficient, __FILE__, __LINE__);
                 
-                Wrappers::Petsc::AXPY(-.5,
+                #ifndef NDEBUG
+                AssertSameNumberingSubset(stiffness, displacement_previous_time_iteration_matrix);
+                #endif // NDEBUG
+                
+                Wrappers::Petsc::AXPY(SAME_NONZERO_PATTERN,
+                                      -.5,
                                       stiffness,
                                       displacement_previous_time_iteration_matrix,
                                       __FILE__, __LINE__);
@@ -201,12 +211,12 @@ namespace HappyHeart
             
             const double factor = 2. / parent::GetTimeManager().GetTimeStep();
 
-            Wrappers::Petsc::AXPY(factor,
+            Wrappers::Petsc::AXPY(SAME_NONZERO_PATTERN, factor,
                                   newest_displacement,
                                   velocity_prev_time_it,
                                   __FILE__, __LINE__);
             
-            Wrappers::Petsc::AXPY(-factor,
+            Wrappers::Petsc::AXPY(SAME_NONZERO_PATTERN, -factor,
                                   former_displacement,
                                   velocity_prev_time_it,
                                   __FILE__, __LINE__);

Sources/ModelInstances/ReactionDiffusion/ReactionDiffusionVariationalFormulation.hxx

@@ -259,7 +259,11 @@ namespace HappyHeart
             auto& system_matrix = this->GetNonCstSystemMatrix(numbering_subset, numbering_subset);
             const auto& capacity_matrix = GetMatrixCapacityPerTimeStep();
             
-            Wrappers::Petsc::AXPY(1., capacity_matrix, system_matrix, __FILE__, __LINE__);
+            #ifndef NDEBUG
+            AssertSameNumberingSubset(capacity_matrix, system_matrix);
+            #endif // NDEBUG
+            
+            Wrappers::Petsc::AXPY(SAME_NONZERO_PATTERN, 1., capacity_matrix, system_matrix, __FILE__, __LINE__);
         }