diff --git a/kratos/future/containers/define_linear_algebra_mpi.h b/kratos/future/containers/define_linear_algebra_mpi.h new file mode 100644 index 000000000000..fc8087bb1b98 --- /dev/null +++ b/kratos/future/containers/define_linear_algebra_mpi.h @@ -0,0 +1,49 @@ +// | / | +// ' / __| _` | __| _ \ __| +// . \ | ( | | ( |\__ ` +// _|\_\_| \__,_|\__|\___/ ____/ +// Multi-Physics +// +// License: BSD License +// Kratos default license: kratos/license.txt +// +// Main authors: Ruben Zorrilla +// + +#pragma once + +// System includes + +// External includes + +// Project includes +#include "containers/distributed_csr_matrix.h" +#include "containers/distributed_sparse_graph.h" +#include "containers/distributed_system_vector.h" + +namespace Kratos::Future +{ + +///@addtogroup KratosCore +///@{ + +///@name Kratos Classes +///@{ + +struct DistributedLinearAlgebraTraits +{ + using DataType = double; + + using IndexType = std::size_t; + + using MatrixType = DistributedCsrMatrix; + + using VectorType = DistributedSystemVector; + + using SparseGraphType = DistributedSparseGraph; +}; + +///@} +///@} addtogroup block + +} // namespace Kratos::Future diff --git a/kratos/future/containers/define_linear_algebra_serial.h b/kratos/future/containers/define_linear_algebra_serial.h new file mode 100644 index 000000000000..836af8f59aee --- /dev/null +++ b/kratos/future/containers/define_linear_algebra_serial.h @@ -0,0 +1,51 @@ +// | / | +// ' / __| _` | __| _ \ __| +// . \ | ( | | ( |\__ ` +// _|\_\_| \__,_|\__|\___/ ____/ +// Multi-Physics +// +// License: BSD License +// Kratos default license: kratos/license.txt +// +// Main authors: Ruben Zorrilla +// + +#pragma once + +// System includes + +// External includes + +// Project includes +#include "containers/csr_matrix.h" +#include "containers/sparse_contiguous_row_graph.h" +#include "containers/system_vector.h" + +namespace Kratos::Future +{ + +///@addtogroup KratosCore +///@{ + +///@name Kratos Classes +///@{ + +struct SerialLinearAlgebraTraits +{ + using DataType = double; + + using IndexType = std::size_t; + + using MatrixType = CsrMatrix; + + using VectorType = SystemVector; + + using DenseMatrixType = DenseMatrix; //TODO: think about this one + + using SparseGraphType = SparseContiguousRowGraph; +}; + +///@} +///@} addtogroup block + +} // namespace Kratos::Future diff --git a/kratos/future/linear_operators/linear_operator.h b/kratos/future/linear_operators/linear_operator.h new file mode 100644 index 000000000000..c82280f9450a --- /dev/null +++ b/kratos/future/linear_operators/linear_operator.h @@ -0,0 +1,227 @@ +// | / | +// ' / __| _` | __| _ \ __| +// . \ | ( | | ( |\__ ` +// _|\_\_| \__,_|\__|\___/ ____/ +// Multi-Physics +// +// License: BSD License +// Kratos default license: kratos/license.txt +// +// Main authors: Ruben Zorrilla +// + +#pragma once + +// System includes + +// External includes + +// Project includes +#include "includes/model_part.h" + +namespace Kratos::Future +{ + +///@addtogroup KratosCore +///@{ + +///@name Kratos Classes +///@{ + +/** + * @brief + * @class LinearOperator + * @ingroup KratosCore + * @brief Auxiliary container to store a linear operator + * @details + * This class implements a generic linear operator interface to be used + * in matrix-free algorithms. Being a generic interface, it can be derived + * to wrap different types of linear operators, including maxtrix-based ones. + * @tparam TLinearAlgebra The struct containing the linear algebra types + */ +template +class LinearOperator +{ + +public: + ///@name Type Definitions + ///@{ + + /// Pointer definition of LinearOperator + KRATOS_CLASS_POINTER_DEFINITION(LinearOperator); + + /// Vector type definition from template parameter + using VectorType = typename TLinearAlgebra::VectorType; + + /// Matrix type definition from template parameter + using MatrixType = typename TLinearAlgebra::MatrixType; + + /// Data type stored in the system vector + using DataType = typename VectorType::DataType; + + /// Index type used in the system vector + using IndexType = typename VectorType::IndexType; + + ///@} + ///@name Life Cycle + ///@{ + + /** + * @brief Default constructor. + * @details Creates an empty LinearOperator with uninitialized function objects. + */ + LinearOperator() = default; + + /** + * @brief Constructor from parameters. + * @param ThisParameters Parameters containing the linear operator settings + */ + LinearOperator(Parameters ThisParameters) + { + mNumRows = ThisParameters["num_rows"].GetInt(); + mNumCols = ThisParameters["num_cols"].GetInt(); + } + + /// Deleted copy constructor (non-copyable) + LinearOperator(const LinearOperator& rOther) = delete; + + /// Defaulted move constructor + LinearOperator(LinearOperator&& rOther) = default; + + /// Default destructor + virtual ~LinearOperator() = default; + + ///@} + ///@name Operators + ///@{ + + /// Deleted copy assignment operator (non-copyable) + LinearOperator& operator=(const LinearOperator& rOther) = delete; + + /// Defaulted move assignment operator + LinearOperator& operator=(LinearOperator&& rOther) = default; + + ///@} + ///@name Operations + ///@{ + + /** + * @brief Performs the matrix-vector product y = A * x. + * @param rX Input vector x + * @param rY Output vector y + */ + virtual void SpMV( + const VectorType& rX, + VectorType& rY) const + { + KRATOS_ERROR << "SpMV() is not implemented in base LinearOperator class." << std::endl; + } + + /** + * @brief Performs the transposed matrix-vector product y = A^T * x. + * @param rX Input vector x + * @param rY Output vector y + */ + virtual void TransposeSpMV( + const VectorType& rX, + VectorType& rY) const + { + KRATOS_ERROR << "TransposeSpMV() is not implemented in base LinearOperator class." << std::endl; + } + + /** + * @brief Clear the operator data. + * @details Resets the sizes and function objects to null. + */ + virtual void Clear() + { + mNumRows = 0; + mNumCols = 0; + } + + ///@} + ///@name Access + ///@{ + + virtual MatrixType& GetMatrix() + { + KRATOS_ERROR << "GetMatrix() is not implemented in base LinearOperator class." << std::endl; + } + + virtual const MatrixType& GetMatrix() const + { + KRATOS_ERROR << "GetMatrix() is not implemented in base LinearOperator class." << std::endl; + } + + /** + * @brief Set the number of rows. + * @param NumRows Number of rows + */ + void SetNumRows(std::size_t NumRows) + { + mNumRows = NumRows; + } + + /** + * @brief Set the number of columns. + * @param NumCols Number of columns + */ + void SetNumCols(std::size_t NumCols) + { + mNumCols = NumCols; + } + + ///@} + ///@name Inquiry + ///@{ + + /** + * @brief Get the number of rows. + * @return Number of rows of the operator + */ + std::size_t NumRows() const + { + return mNumRows; + } + + /** + * @brief Get the number of columns. + * @return Number of columns of the operator + */ + std::size_t NumCols() const + { + return mNumCols; + } + + /** + * @brief Check if the operator is matrix-free. + * @return True if the operator is matrix-free, false otherwise + */ + virtual bool IsMatrixFree() const + { + return true; + } + + ///@} +private: + + ///@name Member Variables + ///@{ + + /// Number of rows of the operator + std::size_t mNumRows = 0; + + /// Number of columns of the operator + std::size_t mNumCols = 0; + + ///@} +}; // class LinearOperator + +///@} +///@name Input and output +///@{ + +///@} +///@} addtogroup block + +} // namespace Kratos::Future diff --git a/kratos/future/linear_operators/sparse_matrix_linear_operator.h b/kratos/future/linear_operators/sparse_matrix_linear_operator.h new file mode 100644 index 000000000000..55343877d531 --- /dev/null +++ b/kratos/future/linear_operators/sparse_matrix_linear_operator.h @@ -0,0 +1,164 @@ +// | / | +// ' / __| _` | __| _ \ __| +// . \ | ( | | ( |\__ ` +// _|\_\_| \__,_|\__|\___/ ____/ +// Multi-Physics +// +// License: BSD License +// Kratos default license: kratos/license.txt +// +// Main authors: Ruben Zorrilla +// + +#pragma once + +// System includes + +// External includes + +// Project includes +#include "includes/model_part.h" +#include "future/linear_operators/linear_operator.h" + +namespace Kratos::Future +{ + +///@addtogroup KratosCore +///@{ + +///@name Kratos Classes +///@{ + +/** + * @brief + * @class SparseMatrixLinearOperator + * @ingroup KratosCore + * @brief Linear operator handling a sparse matrix + * @details + * This class implements a LinearOperator that handles a sparse matrix, + * allowing to use it in matrix-free algorithms while still storing + * the matrix entries explicitly. The linear operator can be also used + * in matrix-based algorithms by accessing the underlying sparse matrix. + * @tparam TLinearAlgebra The struct containing the linear algebra types + */ +template +class SparseMatrixLinearOperator final : public LinearOperator +{ + +public: + ///@name Type Definitions + ///@{ + + /// Pointer definition of SparseMatrixLinearOperator + KRATOS_CLASS_POINTER_DEFINITION(SparseMatrixLinearOperator); + + /// Vector type definition from template parameter + using VectorType = typename TLinearAlgebra::VectorType; + + /// Matrix type definition from template parameter + using MatrixType = typename TLinearAlgebra::MatrixType; + + ///@} + ///@name Life Cycle + ///@{ + + /** + * @brief Default constructor. + * @details Creates an empty SparseMatrixLinearOperator with uninitialized function objects. + */ + SparseMatrixLinearOperator() = default; + + /** + * @brief Constructor from a CSR matrix. + * Constructs a SparseMatrixLinearOperator that wraps the provided CSR matrix. + * Note that this constructor is only enabled when a CSR matrix type is provided. + * @param rA Reference to the CSR matrix + */ + SparseMatrixLinearOperator(MatrixType& rA) + : LinearOperator() + , mrCsrMatrix(rA) + { + this->SetNumRows(rA.size1()); + this->SetNumCols(rA.size2()); + } + + /// Deleted copy constructor (non-copyable) + SparseMatrixLinearOperator(const SparseMatrixLinearOperator& rOther) = delete; + + /// Defaulted move constructor + SparseMatrixLinearOperator(SparseMatrixLinearOperator&& rOther) = default; + + /// Default destructor + ~SparseMatrixLinearOperator() = default; + + ///@} + ///@name Operators + ///@{ + + /// Deleted copy assignment operator (non-copyable) + SparseMatrixLinearOperator& operator=(const SparseMatrixLinearOperator& rOther) = delete; + + /// Defaulted move assignment operator + SparseMatrixLinearOperator& operator=(SparseMatrixLinearOperator&& rOther) = default; + + ///@} + ///@name Operations + ///@{ + + void SpMV( + const VectorType& rX, + VectorType& rY) const override + { + mrCsrMatrix.SpMV(rX, rY); + } + + void TransposeSpMV( + const VectorType& rX, + VectorType& rY) const override + { + mrCsrMatrix.TransposeSpMV(rX, rY); + } + + ///@} + ///@name Access + ///@{ + + MatrixType& GetMatrix() override + { + return mrCsrMatrix; + } + + const MatrixType& GetMatrix() const override + { + return mrCsrMatrix; + } + + ///@} + ///@name Inquiry + ///@{ + + bool IsMatrixFree() const override + { + return false; + } + + ///@} + private: + + ///@name Member Variables + ///@{ + + /// Reference to the CSR matrix + MatrixType& mrCsrMatrix; + + ///@} +}; // class SparseMatrixLinearOperator + +///@} +///@name Input and output +///@{ + +///@} +///@} addtogroup block + +} // namespace Kratos::Future diff --git a/kratos/future/python/add_linear_operators_to_python.cpp b/kratos/future/python/add_linear_operators_to_python.cpp new file mode 100644 index 000000000000..0be75402b1a6 --- /dev/null +++ b/kratos/future/python/add_linear_operators_to_python.cpp @@ -0,0 +1,54 @@ +// | / | +// ' / __| _` | __| _ \ __| +// . \ | ( | | ( |\__ ` +// _|\_\_| \__,_|\__|\___/ ____/ +// Multi-Physics +// +// License: BSD License +// Kratos default license: kratos/license.txt +// +// Main authors: Ruben Zorrilla +// + +// System includes + +// External includes + +// Project includes +#include "includes/define_python.h" +#include "includes/kratos_parameters.h" + +// Future Extensions +#include "future/containers/define_linear_algebra_serial.h" +#include "future/linear_operators/linear_operator.h" +#include "future/linear_operators/sparse_matrix_linear_operator.h" +#include "future/python/add_linear_operators_to_python.h" + +namespace Kratos::Future::Python +{ + +namespace py = pybind11; + +void AddLinearOperatorsToPython(py::module& m) +{ + using LinearOperatorType = Future::LinearOperator; + py::class_(m, "LinearOperator") + .def(py::init()) + .def("SpMV", &LinearOperatorType::SpMV) + .def("TransposeSpMV", &LinearOperatorType::TransposeSpMV) + .def("Clear", &LinearOperatorType::Clear) + .def("SetNumRows", &LinearOperatorType::SetNumRows) + .def("SetNumCols", &LinearOperatorType::SetNumCols) + .def("GetMatrix", py::overload_cast<>(&LinearOperatorType::GetMatrix), py::return_value_policy::reference) + .def("NumRows", &LinearOperatorType::NumRows) + .def("NumCols", &LinearOperatorType::NumCols) + .def("IsMatrixFree", &LinearOperatorType::IsMatrixFree) + ; + + using SparseMatrixLinearOperatorType = Future::SparseMatrixLinearOperator; + py::class_(m, "SparseMatrixLinearOperator") + .def(py::init()) + ; +} + +} // namespace Kratos::Future::Python diff --git a/kratos/future/python/add_linear_operators_to_python.h b/kratos/future/python/add_linear_operators_to_python.h new file mode 100644 index 000000000000..ef65b2943fff --- /dev/null +++ b/kratos/future/python/add_linear_operators_to_python.h @@ -0,0 +1,27 @@ +// | / | +// ' / __| _` | __| _ \ __| +// . \ | ( | | ( |\__ ` +// _|\_\_| \__,_|\__|\___/ ____/ +// Multi-Physics +// +// License: BSD License +// Kratos default license: kratos/license.txt +// +// Main authors: Ruben Zorrilla +// + +#pragma once + +// System includes + +// External includes +#include + +// Project includes + +namespace Kratos::Future::Python +{ + +void AddLinearOperatorsToPython(pybind11::module& m); + +} // namespace Kratos::Future::Python diff --git a/kratos/future/python/kratos_python.cpp b/kratos/future/python/kratos_python.cpp index adb94cd33ef5..3a9538805a18 100644 --- a/kratos/future/python/kratos_python.cpp +++ b/kratos/future/python/kratos_python.cpp @@ -22,6 +22,7 @@ // Future Extensions #include "future/python/add_containers_to_python.h" #include "future/python/add_io_to_python.h" +#include "future/python/add_linear_operators_to_python.h" #include "future/python/add_linear_solvers_to_python.h" #include "future/python/add_processes_to_python.h" #include "future/python/add_strategies_to_python.h" @@ -37,6 +38,8 @@ void AddFutureToPython(py::module& m) AddIOToPython(m); + AddLinearOperatorsToPython(m); + AddLinearSolversToPython(m); AddProcessesToPython(m); diff --git a/kratos/future/tests/cpp_tests/test_sparse_matrix_linear_operator.cpp b/kratos/future/tests/cpp_tests/test_sparse_matrix_linear_operator.cpp new file mode 100644 index 000000000000..89b26b676d16 --- /dev/null +++ b/kratos/future/tests/cpp_tests/test_sparse_matrix_linear_operator.cpp @@ -0,0 +1,66 @@ +// | / | +// ' / __| _` | __| _ \ __| +// . \ | ( | | ( |\__ ` +// _|\_\_| \__,_|\__|\___/ ____/ +// Multi-Physics +// +// License: BSD License +// Kratos default license: kratos/license.txt +// +// Main authors: Ruben Zorrilla +// + +// System includes + +// External includes + +// Project includes +#include "future/containers/define_linear_algebra_serial.h" +#include "future/linear_operators/sparse_matrix_linear_operator.h" +#include "testing/testing.h" + +namespace Kratos::Testing +{ + +KRATOS_TEST_CASE_IN_SUITE(SparseMatrixLinearOperator, KratosCoreFutureSuite) +{ + // Set the input and output vectors + DenseVector aux_data(5); + aux_data[0] = 3.0; + aux_data[1] = 7.0; + aux_data[2] = 2.0; + aux_data[3] = 5.0; + aux_data[4] = 1.0; + SystemVector<> input_vector(aux_data); + SystemVector<> output_vector(5); + output_vector.SetValue(0.0); + + // Set up the CSR matrix + typename CsrMatrix::MatrixMapType matrix_map; + matrix_map[{0, 0}] = 10; matrix_map[{0, 3}] = -2; + matrix_map[{1, 1}] = 8; matrix_map[{1, 2}] = -1; + matrix_map[{2, 2}] = 5; matrix_map[{2, 1}] = -1; matrix_map[{2, 4}] = -2; + matrix_map[{3, 3}] = 7; matrix_map[{3, 0}] = -2; matrix_map[{3, 4}] = -1; + matrix_map[{4, 4}] = 6; matrix_map[{4, 2}] = -2; matrix_map[{4, 3}] = -1; + CsrMatrix csr_matrix(matrix_map); + + // Set up the linear operator from the CSR matrix + const Future::SparseMatrixLinearOperator linear_operator(csr_matrix); + + // Apply the linear operator to an input vector + linear_operator.SpMV(input_vector, output_vector); + + // Check linear operator features + KRATOS_EXPECT_EQ(linear_operator.NumRows(), 5); + KRATOS_EXPECT_EQ(linear_operator.NumCols(), 5); + KRATOS_EXPECT_FALSE(linear_operator.IsMatrixFree()); + + // Check the output values + KRATOS_EXPECT_NEAR(output_vector[0], 20.0, 1e-12); + KRATOS_EXPECT_NEAR(output_vector[1], 54.0, 1e-12); + KRATOS_EXPECT_NEAR(output_vector[2], 1.0, 1e-12); + KRATOS_EXPECT_NEAR(output_vector[3], 28.0, 1e-12); + KRATOS_EXPECT_NEAR(output_vector[4], -3.0, 1e-12); +} + +} \ No newline at end of file diff --git a/kratos/future/tests/test_KratosFutureCore.py b/kratos/future/tests/test_KratosFutureCore.py index aca372dde639..1b4e6aa08a2a 100644 --- a/kratos/future/tests/test_KratosFutureCore.py +++ b/kratos/future/tests/test_KratosFutureCore.py @@ -6,6 +6,7 @@ # Import the tests or test_classes to create the suites import test_vtu_output +import test_sparse_matrix_linear_operator def AssembleTestSuites(): ''' Populates the test suites to run. @@ -24,6 +25,7 @@ def AssembleTestSuites(): # Create a test suite with the selected tests (Small tests): smallSuite = suites['small'] + smallSuite.addTests(KratosUnittest.TestLoader().loadTestsFromTestCases([test_sparse_matrix_linear_operator.TestSparseMatrixLinearOperator])) smallSuite.addTests(KratosUnittest.TestLoader().loadTestsFromTestCases([test_vtu_output.TestVtuOutput2D])) smallSuite.addTests(KratosUnittest.TestLoader().loadTestsFromTestCases([test_vtu_output.TestVtuOutput3D])) smallSuite.addTests(KratosUnittest.TestLoader().loadTestsFromTestCases([test_vtu_output.TestVtuOutput])) diff --git a/kratos/future/tests/test_sparse_matrix_linear_operator.py b/kratos/future/tests/test_sparse_matrix_linear_operator.py new file mode 100644 index 000000000000..8c50cf893a08 --- /dev/null +++ b/kratos/future/tests/test_sparse_matrix_linear_operator.py @@ -0,0 +1,61 @@ +import KratosMultiphysics +import KratosMultiphysics.KratosUnittest as KratosUnittest + +class TestSparseMatrixLinearOperator(KratosUnittest.TestCase): + def test_SparseMatrixLinearOperator(self): + # Create a simple 2x2 CsrMatrix + # [ 2.0 3.0 ] + # [ 1.0 2.0 ] + G = KratosMultiphysics.SparseContiguousRowGraph(2) + G.AddEntry(0,0) + G.AddEntry(0,1) + G.AddEntry(1,0) + G.AddEntry(1,1) + G.Finalize() + + A = KratosMultiphysics.CsrMatrix(G) + A.BeginAssemble() + A.AssembleEntry(2.0,0,0) + A.AssembleEntry(3.0,0,1) + A.AssembleEntry(1.0,1,0) + A.AssembleEntry(2.0,1,1) + A.FinalizeAssemble() + + # Create the CSR matrix linear operator + lin_op = KratosMultiphysics.Future.SparseMatrixLinearOperator(A) + + # Check operator properties + self.assertEqual(lin_op.NumRows(), 2) + self.assertEqual(lin_op.NumCols(), 2) + self.assertFalse(lin_op.IsMatrixFree()) + + # Test GetMatrix + # Note: GetMatrix returns a reference to the CsrMatrix. + A_ref = lin_op.GetMatrix() + self.assertEqual(A_ref.size1(), 2) + + # Set input vector + x = KratosMultiphysics.SystemVector(2) + x.SetValue(1.0) + + # Test SpMV + y = KratosMultiphysics.SystemVector(2) + y_data = y.Data() + lin_op.SpMV(x, y) + self.assertAlmostEqual(y_data[0], 5.0) + self.assertAlmostEqual(y_data[1], 3.0) + + # Test TransposeSpMV + y_trans = KratosMultiphysics.SystemVector(2) + y_trans_data = y_trans.Data() + lin_op.TransposeSpMV(x, y_trans) + self.assertAlmostEqual(y_trans_data[0], 3.0) + self.assertAlmostEqual(y_trans_data[1], 5.0) + + # Test Clear + lin_op.Clear() + self.assertEqual(lin_op.NumRows(), 0) + self.assertEqual(lin_op.NumCols(), 0) + +if __name__ == '__main__': + KratosUnittest.main()