Merge pull request KratosMultiphysics#733 from KratosMultiphysics/fea…

…ture-fsi-mpi

FSI application upgrade

applications/FSIapplication/custom_python/add_convergence_accelerators_to_python.cpp

@@ -47,6 +47,7 @@ void  AddConvergenceAcceleratorsToPython()
 
     // Aitken convergence accelerator
     class_< AitkenConvergenceAccelerator <TSpace>, bases <BaseConvergenceAcceleratorType>, boost::noncopyable > ("AitkenConvergenceAccelerator", init< double >())
+            .def(init< Parameters& >())
             .def("InitializeSolutionStep", &AitkenConvergenceAccelerator<TSpace>::InitializeSolutionStep)
             .def("UpdateSolution", &AitkenConvergenceAccelerator<TSpace>::UpdateSolution)
             .def("FinalizeNonLinearIteration", &AitkenConvergenceAccelerator<TSpace>::FinalizeNonLinearIteration)
@@ -55,6 +56,7 @@ void  AddConvergenceAcceleratorsToPython()
 
     // MVQN convergence accelerator
     class_< MVQNFullJacobianConvergenceAccelerator <TSpace>, bases <BaseConvergenceAcceleratorType>, boost::noncopyable > ("MVQNFullJacobianConvergenceAccelerator", init< double >())
+            .def(init< Parameters& >())
             .def("InitializeSolutionStep", &MVQNFullJacobianConvergenceAccelerator<TSpace>::InitializeSolutionStep)
             .def("UpdateSolution", &MVQNFullJacobianConvergenceAccelerator<TSpace>::UpdateSolution)
             .def("FinalizeNonLinearIteration", &MVQNFullJacobianConvergenceAccelerator<TSpace>::FinalizeNonLinearIteration)
@@ -63,6 +65,7 @@ void  AddConvergenceAcceleratorsToPython()
 
     // MVQN recursive convergence accelerator
     class_< MVQNRecursiveJacobianConvergenceAccelerator <TSpace>, bases <BaseConvergenceAcceleratorType>, boost::noncopyable > ("MVQNRecursiveJacobianConvergenceAccelerator", init< double, int >())
+            .def(init< Parameters& >())
             .def("Initialize", &MVQNRecursiveJacobianConvergenceAccelerator<TSpace>::Initialize)
             .def("InitializeSolutionStep", &MVQNRecursiveJacobianConvergenceAccelerator<TSpace>::InitializeSolutionStep)
             .def("UpdateSolution", &MVQNRecursiveJacobianConvergenceAccelerator<TSpace>::UpdateSolution)

applications/FSIapplication/custom_python/add_custom_utilities_to_python.cpp

@@ -21,6 +21,7 @@
 #include "custom_utilities/FSI_utils.h"
 #include "custom_utilities/aitken_utils.h"
 #include "custom_utilities/partitioned_fsi_utilities.hpp"
+#include "custom_utilities/nodal_update_utilities.h"
 #include "custom_utilities/variable_redistribution_utility.h"
 
 namespace Kratos
@@ -33,35 +34,63 @@ void  AddCustomUtilitiesToPython()
 {
     using namespace boost::python;
     typedef UblasSpace<double, Matrix, Vector > TSpace;
+    typedef NodalUpdateBaseClass< 2 > NodalUpdateBaseClass2DType;
+    typedef NodalUpdateBaseClass< 3 > NodalUpdateBaseClass3DType;
 
     class_<FSIUtils>("FSIUtils", init<>())
 //		.def("FSIUtils",&FSIUtils::GenerateCouplingElements)
-    .def("CheckPressureConvergence",&FSIUtils::CheckPressureConvergence)
-    .def("StructuralPressurePrediction",&FSIUtils::StructuralPressurePrediction)
-    ;
+        .def("CheckPressureConvergence",&FSIUtils::CheckPressureConvergence)
+        .def("StructuralPressurePrediction",&FSIUtils::StructuralPressurePrediction)
+        ;
 
     class_<AitkenUtils>("AitkenUtils", init<>())
-    .def("ComputeAitkenFactor",&AitkenUtils::ComputeAitkenFactor)
-    .def("ComputeRelaxedDisplacement",&AitkenUtils::ComputeRelaxedDisplacement)
-    ;
+        .def("ComputeAitkenFactor",&AitkenUtils::ComputeAitkenFactor)
+        .def("ComputeRelaxedDisplacement",&AitkenUtils::ComputeRelaxedDisplacement)
+        ;
 
     class_<PartitionedFSIUtilities <TSpace,2>, boost::noncopyable >("PartitionedFSIUtilities2D", init< >())
-    .def("GetInterfaceArea",&PartitionedFSIUtilities<TSpace,2>::GetInterfaceArea)
-    .def("GetInterfaceResidualSize",&PartitionedFSIUtilities<TSpace,2>::GetInterfaceResidualSize)
-    .def("SetInterfaceVectorVariable",&PartitionedFSIUtilities<TSpace,2>::SetInterfaceVectorVariable)
-    .def("SetAndFixInterfaceVectorVariable",&PartitionedFSIUtilities<TSpace,2>::SetAndFixInterfaceVectorVariable)
-    .def("ComputeInterfaceVectorResidual",&PartitionedFSIUtilities<TSpace,2>::ComputeInterfaceVectorResidual)
-    .def("ComputeFluidInterfaceMeshVelocityResidualNorm",&PartitionedFSIUtilities<TSpace,2>::ComputeFluidInterfaceMeshVelocityResidualNorm)
-    ;
+        .def("GetInterfaceArea",&PartitionedFSIUtilities<TSpace,2>::GetInterfaceArea)
+        .def("GetInterfaceResidualSize",&PartitionedFSIUtilities<TSpace,2>::GetInterfaceResidualSize)
+        .def("UpdateInterfaceValues",&PartitionedFSIUtilities<TSpace,2>::UpdateInterfaceValues)
+        .def("ComputeInterfaceVectorResidual",&PartitionedFSIUtilities<TSpace,2>::ComputeInterfaceVectorResidual)
+        .def("ComputeFluidInterfaceMeshVelocityResidualNorm",&PartitionedFSIUtilities<TSpace,2>::ComputeFluidInterfaceMeshVelocityResidualNorm)
+        .def("ComputeAndPrintFluidInterfaceNorms",&PartitionedFSIUtilities<TSpace,2>::ComputeAndPrintFluidInterfaceNorms)
+        .def("ComputeAndPrintStructureInterfaceNorms",&PartitionedFSIUtilities<TSpace,2>::ComputeAndPrintStructureInterfaceNorms)
+        .def("CheckCurrentCoordinatesFluid",&PartitionedFSIUtilities<TSpace,2>::CheckCurrentCoordinatesFluid)
+        .def("CheckCurrentCoordinatesStructure",&PartitionedFSIUtilities<TSpace,2>::CheckCurrentCoordinatesStructure)
+        ;
 
     class_<PartitionedFSIUtilities <TSpace,3>, boost::noncopyable >("PartitionedFSIUtilities3D", init< >())
-    .def("GetInterfaceArea",&PartitionedFSIUtilities<TSpace,3>::GetInterfaceArea)
-    .def("GetInterfaceResidualSize",&PartitionedFSIUtilities<TSpace,3>::GetInterfaceResidualSize)
-    .def("SetInterfaceVectorVariable",&PartitionedFSIUtilities<TSpace,3>::SetInterfaceVectorVariable)
-    .def("SetAndFixInterfaceVectorVariable",&PartitionedFSIUtilities<TSpace,3>::SetAndFixInterfaceVectorVariable)
-    .def("ComputeInterfaceVectorResidual",&PartitionedFSIUtilities<TSpace,3>::ComputeInterfaceVectorResidual)
-    .def("ComputeFluidInterfaceMeshVelocityResidualNorm",&PartitionedFSIUtilities<TSpace,3>::ComputeFluidInterfaceMeshVelocityResidualNorm)
-    ;
+        .def("GetInterfaceArea",&PartitionedFSIUtilities<TSpace,3>::GetInterfaceArea)
+        .def("GetInterfaceResidualSize",&PartitionedFSIUtilities<TSpace,3>::GetInterfaceResidualSize)
+        .def("UpdateInterfaceValues",&PartitionedFSIUtilities<TSpace,3>::UpdateInterfaceValues)
+        .def("ComputeInterfaceVectorResidual",&PartitionedFSIUtilities<TSpace,3>::ComputeInterfaceVectorResidual)
+        .def("ComputeFluidInterfaceMeshVelocityResidualNorm",&PartitionedFSIUtilities<TSpace,3>::ComputeFluidInterfaceMeshVelocityResidualNorm)
+        .def("ComputeAndPrintFluidInterfaceNorms",&PartitionedFSIUtilities<TSpace,3>::ComputeAndPrintFluidInterfaceNorms)
+        .def("ComputeAndPrintStructureInterfaceNorms",&PartitionedFSIUtilities<TSpace,3>::ComputeAndPrintStructureInterfaceNorms)
+        .def("CheckCurrentCoordinatesFluid",&PartitionedFSIUtilities<TSpace,3>::CheckCurrentCoordinatesFluid)
+        .def("CheckCurrentCoordinatesStructure",&PartitionedFSIUtilities<TSpace,3>::CheckCurrentCoordinatesStructure)
+        ;
+
+    class_< NodalUpdateBaseClass <2>, boost::noncopyable > ("BaseNodalUpdate2D", init< >())
+        .def("UpdateMeshTimeDerivatives", &NodalUpdateBaseClass<2>::UpdateMeshTimeDerivatives)
+        .def("UpdateTimeDerivativesOnInterface", &NodalUpdateBaseClass<2>::UpdateTimeDerivativesOnInterface)
+        ;
+
+    class_< NodalUpdateBaseClass <3>, boost::noncopyable > ("BaseNodalUpdate3D", init< >())
+        .def("UpdateMeshTimeDerivatives", &NodalUpdateBaseClass<3>::UpdateMeshTimeDerivatives)
+        .def("UpdateTimeDerivativesOnInterface", &NodalUpdateBaseClass<3>::UpdateTimeDerivativesOnInterface)
+        ;
+
+    class_< NodalUpdateNewmark <2>, bases <NodalUpdateBaseClass2DType>, boost::noncopyable > ("NodalUpdateNewmark2D", init< const double >())
+        .def("UpdateMeshTimeDerivatives", &NodalUpdateNewmark<2>::UpdateMeshTimeDerivatives)
+        .def("UpdateTimeDerivativesOnInterface", &NodalUpdateNewmark<2>::UpdateTimeDerivativesOnInterface)
+        ;
+
+    class_< NodalUpdateNewmark <3>, bases <NodalUpdateBaseClass3DType>, boost::noncopyable > ("NodalUpdateNewmark3D", init< const double >())
+        .def("UpdateMeshTimeDerivatives", &NodalUpdateNewmark<3>::UpdateMeshTimeDerivatives)
+        .def("UpdateTimeDerivativesOnInterface", &NodalUpdateNewmark<3>::UpdateTimeDerivativesOnInterface)
+        ;
 
     typedef void (*DistributePointDoubleType)(ModelPart&, const Variable< double >&, const Variable< double >&, double, unsigned int);
     typedef void (*DistributePointArrayType)(ModelPart&, const Variable< array_1d<double,3> >&, const Variable< array_1d<double,3> >&,double, unsigned int);

applications/FSIapplication/custom_python/fsi_python_application.cpp

@@ -43,25 +43,6 @@ BOOST_PYTHON_MODULE(KratosFSIApplication)
     AddMappersToPython();
     AddConvergenceAcceleratorsToPython();
 
-    //registering variables in python
-    KRATOS_REGISTER_IN_PYTHON_VARIABLE(CONVERGENCE_ACCELERATOR_ITERATION);
-    KRATOS_REGISTER_IN_PYTHON_VARIABLE(MAPPER_SCALAR_PROJECTION_RHS);
-    KRATOS_REGISTER_IN_PYTHON_VARIABLE(FICTITIOUS_FLUID_DENSITY);
-    KRATOS_REGISTER_IN_PYTHON_VARIABLE(SCALAR_PROJECTED);
-    KRATOS_REGISTER_IN_PYTHON_VARIABLE(FSI_INTERFACE_RESIDUAL_NORM);
-    KRATOS_REGISTER_IN_PYTHON_VARIABLE(FSI_INTERFACE_MESH_RESIDUAL_NORM);
-//	KRATOS_REGISTER_IN_PYTHON_VARIABLE(PRESSURE_OLD_IT);
-//	KRATOS_REGISTER_IN_PYTHON_VARIABLE(IS_INTERFACE);
-
-    KRATOS_REGISTER_IN_PYTHON_3D_VARIABLE_WITH_COMPONENTS(MAPPER_VECTOR_PROJECTION_RHS);
-    KRATOS_REGISTER_IN_PYTHON_3D_VARIABLE_WITH_COMPONENTS(VAUX_EQ_TRACTION);
-    KRATOS_REGISTER_IN_PYTHON_3D_VARIABLE_WITH_COMPONENTS(VECTOR_PROJECTED);
-    KRATOS_REGISTER_IN_PYTHON_3D_VARIABLE_WITH_COMPONENTS(RELAXED_DISP);
-    KRATOS_REGISTER_IN_PYTHON_3D_VARIABLE_WITH_COMPONENTS(FSI_INTERFACE_RESIDUAL);
-    KRATOS_REGISTER_IN_PYTHON_3D_VARIABLE_WITH_COMPONENTS(FSI_INTERFACE_MESH_RESIDUAL);
-    KRATOS_REGISTER_IN_PYTHON_3D_VARIABLE_WITH_COMPONENTS(POSITIVE_MAPPED_VECTOR_VARIABLE);
-    KRATOS_REGISTER_IN_PYTHON_3D_VARIABLE_WITH_COMPONENTS(NEGATIVE_MAPPED_VECTOR_VARIABLE);
-
 }
 
 

applications/FSIapplication/custom_utilities/aitken_convergence_accelerator.hpp

@@ -14,12 +14,12 @@
 /* System includes */
 
 /* External includes */
-#include "boost/smart_ptr.hpp"
 #include "utilities/math_utils.h"
 
 /* Project includes */
 #include "includes/define.h"
 #include "includes/variables.h"
+#include "includes/kratos_parameters.h"
 #include "convergence_accelerator.hpp"
 
 namespace Kratos
@@ -69,6 +69,22 @@ class AitkenConvergenceAccelerator: public ConvergenceAccelerator<TSpace>
      * Constructor.
      * Aitken convergence accelerator
      */
+    AitkenConvergenceAccelerator(Parameters& rConvAcceleratorParameters)
+    {
+        Parameters aitken_default_parameters(R"(
+        {
+            "solver_type"       : "Relaxation",
+            "acceleration_type" : "Aitken",
+            "w_0"               : 0.825  
+        }
+        )");
+
+        rConvAcceleratorParameters.ValidateAndAssignDefaults(aitken_default_parameters);
+
+        mOmega_0 = rConvAcceleratorParameters["w_0"].GetDouble();
+    }
+
+
     AitkenConvergenceAccelerator(double rOmegaInitial = 0.825)
     {
         mOmega_0 = rOmegaInitial;
@@ -128,28 +144,23 @@ class AitkenConvergenceAccelerator: public ConvergenceAccelerator<TSpace>
     {
         KRATOS_TRY;
 
-        // (*mpResidualVector_1) = rResidualVector;
         VectorPointerType pAux(new VectorType(rResidualVector));
         std::swap(mpResidualVector_1, pAux);
 
         if (mConvergenceAcceleratorIteration == 1)
         {
-            // rIterationGuess += mOmega_0*mResidualVector_1;
             TSpace::UnaliasedAdd(rIterationGuess, mOmega_0, *mpResidualVector_1);
         }
         else
         {
             VectorType Aux1minus0(*mpResidualVector_1);                  // Auxiliar copy of mResidualVector_1
             TSpace::UnaliasedAdd(Aux1minus0, -1.0, *mpResidualVector_0); // mResidualVector_1 - mResidualVector_0
 
-            // double num = inner_prod(mResidualVector_0,(mResidualVector_1-mResidualVector_0));
-            // double den = inner_prod((mResidualVector_1-mResidualVector_0),(mResidualVector_1-mResidualVector_0));
             double den = TSpace::Dot(Aux1minus0, Aux1minus0);
             double num = TSpace::Dot(*mpResidualVector_0, Aux1minus0);
 
             mOmega_1 = -mOmega_0*(num/den);
 
-            // rIterationGuess += mOmega_1*mResidualVector_1;
             TSpace::UnaliasedAdd(rIterationGuess, mOmega_1, *mpResidualVector_1);
             mOmega_0 = mOmega_1;
         }
@@ -164,8 +175,8 @@ class AitkenConvergenceAccelerator: public ConvergenceAccelerator<TSpace>
     {
         KRATOS_TRY;
 
-        // *mResidualVector_0 = *mResidualVector_1;
-        mpResidualVector_0 = mpResidualVector_1;
+        // mpResidualVector_0 = mpResidualVector_1;
+        std::swap(mpResidualVector_0, mpResidualVector_1);
         mConvergenceAcceleratorIteration += 1;
 
         KRATOS_CATCH( "" );

applications/FSIapplication/custom_utilities/mvqn_convergence_accelerator.hpp

@@ -14,13 +14,13 @@
 /* System includes */
 
 /* External includes */
-#include "boost/smart_ptr.hpp"
 #include "utilities/math_utils.h"
 
 /* Project includes */
 #include "includes/define.h"
 #include "includes/variables.h"
 #include "includes/ublas_interface.h"
+#include "includes/kratos_parameters.h"
 #include "convergence_accelerator.hpp"
 
 namespace Kratos
@@ -73,6 +73,22 @@ class MVQNFullJacobianConvergenceAccelerator: public ConvergenceAccelerator<TSpa
      * Constructor.
      * Aitken convergence accelerator
      */
+    MVQNFullJacobianConvergenceAccelerator( Parameters &rConvAcceleratorParameters )
+    {
+        Parameters mvqn_default_parameters(R"(
+        {
+            "solver_type" : "MVQN",
+            "w_0"         : 0.825
+        }
+        )");
+
+        rConvAcceleratorParameters.ValidateAndAssignDefaults(mvqn_default_parameters);
+
+        mOmega_0 = rConvAcceleratorParameters["w_0"].GetDouble();
+        mConvergenceAcceleratorIteration = 0;
+        mConvergenceAcceleratorFirstCorrectionPerformed = false;
+    }
+
     MVQNFullJacobianConvergenceAccelerator( double rOmegaInitial = 0.825 )
     {
         mOmega_0 = rOmegaInitial;

applications/FSIapplication/custom_utilities/mvqn_recursive_convergence_accelerator.hpp

@@ -180,43 +180,40 @@ class JacobianEmulator
         {
             if (mpOldJacobianEmulator != nullptr) // If it is available, consider the previous step Jacobian
             {
-                // mpOldJacobianEmulator->ApplyJacobian(rWorkVector, rProjectedVector);
                 mpOldJacobianEmulator->ApplyJacobian(pWorkVector, pProjectedVector);
             }
             else // When the JacobianEmulator has no PreviousJacobianEmulator consider minus the identity matrix as inverse Jacobian
             {
-                // TSpace::Assign(rProjectedVector,-1.0,rWorkVector);
                 TSpace::Assign(*pProjectedVector, -1.0, *pWorkVector);
             }
         }
         else
         {
+            const unsigned int previous_iterations = mJacobianObsMatrixV.size();
+            const unsigned int residual_size = TSpace::Size(mJacobianObsMatrixV[0]);
 
-            VectorPointerType pY(new VectorType(mJacobianObsMatrixV[0]));
-            VectorPointerType pW(new VectorType(mJacobianObsMatrixV[0]));
-            VectorPointerType pzQR(new VectorType(mJacobianObsMatrixV.size()));
-            MatrixPointerType pAuxMatQR(new MatrixType(TSpace::Size(mJacobianObsMatrixV[0]), mJacobianObsMatrixV.size()));
+            VectorPointerType pY(new VectorType(residual_size));
+            VectorPointerType pW(new VectorType(residual_size));
+            VectorPointerType pzQR(new VectorType(previous_iterations));
+            MatrixPointerType pAuxMatQR(new MatrixType(residual_size, previous_iterations));
 
             // Loop to store a std::vector<VectorType> type as Matrix type
-            for (unsigned int i = 0; i < TSpace::Size(mJacobianObsMatrixV[0]); i++)
+            for (unsigned int i = 0; i < residual_size; ++i)
             {
-                for (unsigned int j = 0; j < mJacobianObsMatrixV.size(); j++)
+                for (unsigned int j = 0; j < previous_iterations; ++j)
                 {
-                    // auxMatQR(i,j) = mJacobianObsMatrixV[j](i);
                     (*pAuxMatQR)(i,j) = mJacobianObsMatrixV[j](i);
                 }
             }
 
             VectorPointerType pWorkVectorCopy(new VectorType(*pWorkVector));
 
             // QR decomposition to compute ((V_k.T*V_k)^-1)*V_k.T*r_k
-            // TODO: Implement an if(!frozen) to avoid the recomputation of the inverse matrix each time the OldJacobianEmulator is called
-            mQR_decomposition.compute(TSpace::Size(mJacobianObsMatrixV[0]), mJacobianObsMatrixV.size(), &(*pAuxMatQR)(0,0));
+            mQR_decomposition.compute(residual_size, previous_iterations, &(*pAuxMatQR)(0,0));
             mQR_decomposition.solve(&(*pWorkVectorCopy)(0), &(*pzQR)(0));
 
-            // TODO: PARALLELIZE THIS OPERATION (it cannot be done with TSpace::Dot beacuse the obs matrices are stored by columns)
             TSpace::SetToZero(*pY);
-            for (unsigned int j = 0; j < mJacobianObsMatrixV.size(); j++)
+            for (unsigned int j = 0; j < previous_iterations; ++j)
             {
                 TSpace::UnaliasedAdd(*pY, (*pzQR)(j), mJacobianObsMatrixV[j]);
             }
@@ -236,7 +233,7 @@ class JacobianEmulator
 
             // w = W_k*z
             TSpace::SetToZero(*pW);
-            for (unsigned int j = 0; j < mJacobianObsMatrixV.size(); j++)
+            for (unsigned int j = 0; j < previous_iterations; ++j)
             {
                 TSpace::UnaliasedAdd(*pW, (*pzQR)(j), mJacobianObsMatrixW[j]);
             }
@@ -433,6 +430,25 @@ class MVQNRecursiveJacobianConvergenceAccelerator: public ConvergenceAccelerator
      * Constructor.
      * MVQN convergence accelerator
      */
+    MVQNRecursiveJacobianConvergenceAccelerator( Parameters &rConvAcceleratorParameters )
+    {
+        Parameters mvqn_recursive_default_parameters(R"(
+        {
+            "solver_type" : "MVQN",
+            "w_0"         : 0.825,
+            "buffer_size" : 10
+        }
+        )");
+
+        rConvAcceleratorParameters.ValidateAndAssignDefaults(mvqn_recursive_default_parameters);
+
+        mOmega_0 = rConvAcceleratorParameters["w_0"].GetDouble();
+        mJacobianBufferSize = rConvAcceleratorParameters["buffer_size"].GetInt();
+        mConvergenceAcceleratorStep = 0;
+        mConvergenceAcceleratorIteration = 0;
+        mConvergenceAcceleratorFirstCorrectionPerformed = false;
+    }
+
     MVQNRecursiveJacobianConvergenceAccelerator( double rOmegaInitial = 0.825, unsigned int rJacobianBufferSize = 10 )
     {
         mOmega_0 = rOmegaInitial;