#723 FSI/Newton: solver to solve differential of solid equation introduced.

Data/Lua/demo_input_fsi_ei.lua

@@ -264,7 +264,7 @@ EssentialBoundaryCondition3 = {
 Mesh1 = {
     -- Path of the mesh file to use.
     -- Expected format: "VALUE"
-    mesh = "${HOME}/Codes/HappyHeart/Data/Mesh/ThTotal2D.mesh",
+    mesh = "${HOME}/Codes/HappyHeart/Data/Mesh/ThTotal2Dsmall.mesh",
     
     -- Format of the input mesh.
     -- Expected format: "VALUE"
@@ -622,7 +622,7 @@ Petsc = {
     -- Absolute tolerance
     -- Expected format: {VALUE1, VALUE2, ...}
     -- Constraint: v > 0.
-    absoluteTolerance = {1e-10},
+    absoluteTolerance = 1e-10,
     
     -- gmresStart
     -- Expected format: {VALUE1, VALUE2, ...}
@@ -644,12 +644,12 @@ Petsc = {
     -- Relative tolerance
     -- Expected format: {VALUE1, VALUE2, ...}
     -- Constraint: v > 0.
-    relativeTolerance = {1e-6},
+    relativeTolerance = 1e-6,
     
     -- List of solver: { chebychev cg gmres preonly bicg python };
     -- Expected format: {"VALUE1", "VALUE2", ...}
     -- Constraint: ops_in(v, {'Mumps', 'Umfpack'})
-    solver = {'Umfpack'},
+    solver = 'Umfpack',
     
 } -- Petsc
 

Data/Lua/demo_input_hyperelasticity.lua

@@ -581,7 +581,7 @@ TransientSource2 = {
 	-- Value of the transient source in the case nature is 'constant'(and also initial value if nature is 
 	-- 'at_quadrature_point'; irrelevant otherwise). 
 	-- Expected format: {VALUE1, VALUE2, ...}
-	scalar_value = { 0. , 0.005, 0.},
+	scalar_value = { 0. , 0.5, 0.},
 
 	-- Value of the transient source  in the case nature is 'lua_function'(and also initial value if nature is 
 	-- 'at_quadrature_point'; irrelevant otherwise).

Sources/ModelInstances/FSI_EI/Aitken/ModelForward.hxx

@@ -200,9 +200,8 @@ namespace HappyHeart
                 if (mpi.IsRootProcessor())
                     std::cout << "Numbering = " << numbering_subset.GetUniqueId() << std::endl;
                 
-                implicit_step_fluid_formulation.DefineFormulation();
-                
-                
+                implicit_step_fluid_formulation.DefineFormulation(differential::no);
+                                
                 implicit_step_fluid_formulation.GetSystemRhs(numbering_subset).template Print<MpiScale::processor_wise>(mpi, aitken_output_dir + "implicit_step_fluid_rhs.hhdata", __FILE__, __LINE__);
                 
                 implicit_step_fluid_formulation.template ApplyEssentialBoundaryCondition<OperatorNS::Nature::nonlinear>(numbering_subset,

Sources/ModelInstances/FSI_EI/Aitken/main_aitken.cpp

@@ -58,7 +58,7 @@ int main(int argc, char ** argv)
             FSI_EINS::Policy::Hyperelasticity
             <
                 HyperelasticityNS::LawNS::CiarletGeymonat,
-                HyperelasticityNS::TimeScheme::midpoint
+                HyperelasticityNS::TimeScheme::half_sum
             >;
         
 #else

Sources/ModelInstances/FSI_EI/Newton/Model.hpp

@@ -233,6 +233,9 @@ namespace HappyHeart
             GlobalVector& GetNonCstSolidDisplacementBeforeSoFOnInterface() noexcept;
 
             
+            //! Access to the linear solver used to solve the diffferential of the solid variational formulation.
+            Wrappers::Petsc::Snes& GetNonCstDifferentialSolidSolver() noexcept;
+            
             
         private:
             
@@ -266,6 +269,8 @@ namespace HappyHeart
             GlobalVector::unique_ptr work_vector_2_ = nullptr;
             GlobalVector::unique_ptr work_vector_3_ = nullptr;
             
+            GlobalVector::unique_ptr solid_displacement_inside_dh_ = nullptr;
+            
             /*!
              * \brief Interpolation matrix for solid velocity -> fluid velocity.
              *
@@ -277,6 +282,8 @@ namespace HappyHeart
             ConformInterpolatorNS::SubsetOrSuperset::unique_ptr solid_to_fsi_interpolator_ = nullptr;
             
             ConformInterpolatorNS::P1_to_P2::unique_ptr solid_to_fluid_velocity_on_fsi_interpolator_ = nullptr;
+            
+            Wrappers::Petsc::Snes::unique_ptr differential_solid_solver_ = nullptr;
 
         };
         

Sources/ModelInstances/FSI_EI/Newton/Model.hxx

@@ -119,6 +119,9 @@ namespace HappyHeart
                 auto& solid_displacement_before_SoF = GetNonCstSolidDisplacementBeforeSoF();
                 HappyHeart::AllocateGlobalVector(god_of_dof, solid_displacement_before_SoF);
                 
+                
+                solid_displacement_inside_dh_ = std::make_unique<GlobalVector>(solid_displacement_before_SoF);
+                
                 solid_displacement_evolution_in_SoF_ = std::make_unique<GlobalVector>(solid_displacement_before_SoF);
                 
                 solid_displacement_before_SoF_on_interface_ =
@@ -244,6 +247,7 @@ namespace HappyHeart
                 solid_to_fsi_interpolator_->Init();
             }
             
+            differential_solid_solver_ = InitSolver(mpi, input_parameter_data);            
         }
         
              

Sources/ModelInstances/FSI_EI/Newton/ModelAccessors.hxx

@@ -160,7 +160,14 @@ namespace HappyHeart
         }
 
     
-     
+        template<class SolidVariationalFormulationPolicyT>
+        Wrappers::Petsc::Snes& Model<SolidVariationalFormulationPolicyT>::GetNonCstDifferentialSolidSolver() noexcept
+        {
+            assert(!(!differential_solid_solver_));
+            return *differential_solid_solver_;
+        }
+        
+        
     } // namespace FSI_EINS
     
 

Sources/ModelInstances/FSI_EI/Newton/ModelForward.hxx

@@ -47,8 +47,31 @@ namespace HappyHeart
             
             UpdateFluidDirichletCondition(differential::yes);
             ImplicitFluidStep(differential::yes);
+            auto& fluid_varf = GetNonCstImplicitStepFluidVariationalFormulation();
+            const auto& differential_fluid_residual = fluid_varf.ComputeResidual();
             
-//            const auto& fluid_varf = GetImplicitStepFluidVariationalFormulation();
+            const auto& solid_varf = SolidVariationalFormulationPolicyT::GetVariationalFormulation();
+            const auto& solid_numbering_subset = solid_varf.GetNumberingSubset();
+            const auto& solid_tangent = solid_varf.GetSystemMatrix(solid_numbering_subset, solid_numbering_subset);
+            
+           // auto& solid_residual = *solid_displacement_inside_dh_;
+            
+            auto& solid_residual = GetNonCstSolidResidual();
+            Wrappers::Petsc::MatMultTranspose(GetInterpolationMatrixSolid2FluidVelocity(),
+                                              differential_fluid_residual,
+                                              solid_residual,
+                                              __FILE__, __LINE__);
+            
+//            GetNonCstDifferentialSolidSolver().
+            
+            // SOLVE LINEAR HERE!
+            solid_residual.View(mpi, __FILE__, __LINE__);
+            
+            
+
+            
+            
+//
 //            const auto& sol = fluid_varf.GetSystemSolution(fluid_varf.GetNumberingSubset());
 //            sol.View(mpi, __FILE__, __LINE__);
             

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

@@ -58,11 +58,11 @@ namespace HappyHeart
             void ComputeDynamicTimeSchemeTangentContribution(const VectorsAndMatrices<HyperelasticityNS::TimeScheme::half_sum>& vm,
                                                              GlobalMatrix& system_matrix)
             {
-                Wrappers::Petsc::AXPY(1.,
+                Wrappers::Petsc::AXPY(0.5,
                                       vm.GetMatrixNewStiffness(),
                                       system_matrix, __FILE__, __LINE__);
                 
-                Wrappers::Petsc::AXPY(1.,
+                Wrappers::Petsc::AXPY(0.5,
                                       vm.GetMatrixStiffnessPreviousTimeIteration(),
                                       system_matrix, __FILE__, __LINE__);
             }
@@ -71,12 +71,12 @@ namespace HappyHeart
             void ComputeDynamicTimeSchemeResidualContribution(const VectorsAndMatrices<HyperelasticityNS::TimeScheme::half_sum>& vm,
                                                               GlobalVector& system_rhs)
             {
-                Wrappers::Petsc::AXPY(1.,
+                Wrappers::Petsc::AXPY(0.5,
                                       vm.GetVectorNewStiffness(),
                                       system_rhs,
                                       __FILE__, __LINE__);
 
-                Wrappers::Petsc::AXPY(1.,
+                Wrappers::Petsc::AXPY(0.5,
                                       vm.GetVectorStiffnessPreviousTimeIteration(),
                                       system_rhs,
                                       __FILE__, __LINE__);