Smoothers and Projections seem working

integration_tests/CMakeLists.txt

@@ -47,6 +47,15 @@ target_include_directories(int_test_heat_inhomogenious_CG
         INTERFACE ${CMAKE_CURRENT_SOURCE_DIR})
 target_compile_features(int_test_heat_inhomogenious_CG PRIVATE cxx_std_17)
 
+add_executable(int_test_multigrid int_test_multigrid.cpp
+        heat_solver.cpp
+        heat_solver.hpp)
+target_link_libraries(int_test_multigrid PRIVATE Catch2::Catch2WithMain vtint mfem_helper ${MFEM_LIBRARIES})
+target_include_directories(int_test_multigrid
+        PUBLIC ${MFEM_INCLUDE_DIRS}
+        INTERFACE ${CMAKE_CURRENT_SOURCE_DIR})
+target_compile_features(int_test_multigrid PRIVATE cxx_std_17)
+
 add_executable(int_test_ode int_test_ode.cpp)
 target_link_libraries(int_test_ode PRIVATE Catch2::Catch2WithMain vtint mfem_helper ${MFEM_LIBRARIES})
 target_include_directories(int_test_ode

integration_tests/heat_solver.cpp

@@ -51,11 +51,13 @@ get_inhomogenious_rhs(mfem::FiniteElementSpace &T_test, std::function<mfem::Vect
     return b_I;
 }
 
-double
-integration_tests::solve_heat(const int nx, const int ny, const int nt, const int order_x, const int order_t,
-                              std::function<double(double, mfem::Vector const&)> c1, std::function<double(double, mfem::Vector const&)> c2,
-                              std::function<double(double, mfem::Vector const&)> f_analytical,
-                              std::function<double(double, mfem::Vector const&)> u_analytical) {
+std::tuple<Eigen::SparseMatrix<double>, Eigen::VectorXd, vtint::block_vector_t, mfem::GridFunction>
+integration_tests::assemble_heat_system_2d_DG(const int nx, const int ny, const int nt, const int order_x,
+                                                   const int order_t,
+                                                   std::function<double(double, const mfem::Vector &)> c1,
+                                                   std::function<double(double, const mfem::Vector &)> c2,
+                                                   std::function<double(double, const mfem::Vector &)> f_analytical,
+                                                   std::function<double(double, const mfem::Vector &)> u_analytical) {
     //construct meshes
     mfem::Mesh mesh_x = mfem::Mesh::MakeCartesian2D(nx, ny, mfem::Element::TRIANGLE, true, 1.0, 1.0);
     mfem::Mesh mesh_t = mfem::Mesh::MakeCartesian1D(nt, 1.0);
@@ -150,6 +152,30 @@ integration_tests::solve_heat(const int nx, const int ny, const int nt, const in
     Eigen::VectorXd const b = vtint::homogenize_rhs(A_h, A_rhs, b_I.to_vector(),
                                                     u_initial.get_slice(vtint::range(1, u_initial.get_n_blocks())).to_vector(), ess_dof_q);
 
+    return std::make_tuple(A_h, b, u_initial, u0);
+}
+
+
+double
+integration_tests::solve_heat(const int nx, const int ny, const int nt, const int order_x, const int order_t,
+                              std::function<double(double, mfem::Vector const&)> c1, std::function<double(double, mfem::Vector const&)> c2,
+                              std::function<double(double, mfem::Vector const&)> f_analytical,
+                              std::function<double(double, mfem::Vector const&)> u_analytical) {
+    //construct meshes
+    mfem::Mesh mesh_x = mfem::Mesh::MakeCartesian2D(nx, ny, mfem::Element::TRIANGLE, true, 1.0, 1.0);
+    mfem::Mesh mesh_t = mfem::Mesh::MakeCartesian1D(nt, 1.0);
+
+    //construct finite element spaces
+    mfem::H1_FECollection fec_x(order_x, mesh_x.Dimension());
+    mfem::H1_FECollection fec_t(order_t, mesh_t.Dimension());
+    mfem::L2_FECollection fec_t_test(order_t-1, mesh_t.Dimension());
+    mfem::FiniteElementSpace X(&mesh_x, &fec_x);
+    mfem::FiniteElementSpace T(&mesh_t, &fec_t);
+    mfem::FiniteElementSpace T_test(&mesh_t, &fec_t_test);
+
+    //assemble FE Matrices
+    auto [A_h, b, u_initial, u0] = assemble_heat_system_2d_DG(nx, ny, nt, order_x, order_t, c1, c2, f_analytical, u_analytical);
+
     //solve system with eigen LU solver
     Eigen::SparseLU<Eigen::SparseMatrix<double>> solver;
     solver.compute(A_h);
@@ -464,3 +490,4 @@ double integration_tests::solve_heat_CG(const int nx, const int ny, const int nt
 
     return l2_error;
 }
+

integration_tests/heat_solver.hpp

@@ -9,10 +9,13 @@
 #include <mfem.hpp>
 #include "block_vector_t.hpp"
 
+#include <Eigen/SparseCore>
+
 namespace integration_tests {
     double
     solve_heat(const int nx, const int ny, const int nt, const int order_x, const int order_t,
-               std::function<double(double, mfem::Vector const&)> c1, std::function<double(double, mfem::Vector const&)> c2,
+               std::function<double(double, mfem::Vector const &)> c1,
+               std::function<double(double, mfem::Vector const &)> c2,
                std::function<double(double, mfem::Vector const &)> f_analytical,
                std::function<double(double, mfem::Vector const &)> u_analytical);
 
@@ -29,6 +32,14 @@ namespace integration_tests {
     solve_heat_CG(const int nx, const int ny, const int nt, const int order_x, const int order_t,
                   std::function<double(double, mfem::Vector const &)> f_analytical,
                   std::function<double(double, mfem::Vector const &)> u_analytical);
+
+    std::tuple<Eigen::SparseMatrix<double>, Eigen::VectorXd, vtint::block_vector_t, mfem::GridFunction>
+    assemble_heat_system_2d_DG(const int nx, const int ny, const int nt, const int order_x, const int order_t,
+                               std::function<double(double, mfem::Vector const &)> c1,
+                               std::function<double(double, mfem::Vector const &)> c2,
+                               std::function<double(double, mfem::Vector const &)> f_analytical,
+                               std::function<double(double, mfem::Vector const &)> u_analytical);
+
 } // integration_tests
 
 #endif //VTINT_HEAT_SOLVER_HPP

integration_tests/int_test_multigrid.cpp

+//
+// Created by mreic on 02.09.2024.
+//
+
+#include <catch2/catch_test_macros.hpp>
+#include <catch2/matchers/catch_matchers_floating_point.hpp>
+
+#include "simple_multigrid.hpp"
+#include "heat_solver.hpp"
+#include "homogenization.hpp"
+
+#include <Eigen/SparseLU>
+#include "mfem.hpp"
+
+TEST_CASE("VTINT: Multigrid", "[vtint][multigrid]") {
+    int const n = 2;
+    mfem::Mesh mesh = mfem::Mesh::MakeCartesian1D(n);
+    mfem::H1_FECollection fecX(1, mesh.Dimension());
+    mfem::L2_FECollection fecY(0, mesh.Dimension());
+    mfem::FiniteElementSpace X(&mesh, &fecX);
+    mfem::FiniteElementSpace Y(&mesh, &fecY);
+
+    SECTION("See how projections work"){
+
+        mesh.UniformRefinement();
+        mfem::FiniteElementSpace X1(&mesh, &fecX);
+        mfem::FiniteElementSpace Y1(&mesh, &fecY);
+
+        mfem::OperatorHandle Px;
+        X1.GetTransferOperator(X,Px);
+
+        mfem::SparseMatrix Pm = *Px.As<mfem::SparseMatrix>();
+        Pm.Finalize();
+        Pm.ToDenseMatrix()->Print(std::cout);
+
+        mfem::OperatorHandle Py;
+        Y1.GetTransferOperator(Y,Py);
+
+        mfem::SparseMatrix Pyn = *Py.As<mfem::SparseMatrix>();
+        Pyn.Finalize();
+        Pyn.ToDenseMatrix()->Print(std::cout);
+        REQUIRE(1 == 1);
+    }SECTION("Generate Grid Hierarchy"){
+        vtint::grid_hierarchy_t grid_hierarchy(3, X, Y);
+        REQUIRE(1 == 1);
+    }
+
+
+    REQUIRE(1 == 1);
+}
+
+
+TEST_CASE("VTINT: Multigrid", "[vtint][multigrid][smoother]"){
+    int const n = 2;
+    mfem::Mesh mesh = mfem::Mesh::MakeCartesian1D(n);
+    mfem::H1_FECollection fecX(1, mesh.Dimension());
+    mfem::L2_FECollection fecY(0, mesh.Dimension());
+    mfem::FiniteElementSpace X(&mesh, &fecX);
+    mfem::FiniteElementSpace Y(&mesh, &fecY);
+
+    static int constexpr order_x = 1;
+    static int constexpr order_t = 1;
+    static int constexpr nx = 10;
+    static int constexpr ny = nx;
+    static int constexpr nt = 10;
+
+    auto c1 = [](double t, mfem::Vector const& p){
+        return 1.0;
+    };
+    auto c2 = [](double t, mfem::Vector const& p){
+        return 1.0;
+    };
+    auto u_analytical = [](double t, mfem::Vector const& p){
+        double const& x = p[0];
+        double const& y = p[1];
+        return std::sin(M_PI*x)*std::sin(M_PI*y)*t*t;
+    };
+    auto f_analytical = [](double t, mfem::Vector const& p){
+        double const& x = p[0];
+        double const& y = p[1];
+        return 2.0*std::sin(M_PI*x)*std::sin(M_PI*y)*t+2.0*M_PI*M_PI*std::sin(M_PI*x)*std::sin(M_PI*y)*t*t;
+    };
+
+    //take a simple heat DG system as base
+    auto [A_h, b_h, u_initial, u_at_t0] = integration_tests::assemble_heat_system_2d_DG(nx, ny, nt, order_x, order_t, c1, c2, f_analytical, u_analytical);
+    Eigen::Index block_size = u_initial.get_block(0).size();
+    Eigen::Index rows = A_h.rows();
+    Eigen::Index cols = A_h.cols();
+    REQUIRE(rows == cols);
+    REQUIRE(rows % block_size == 0);
+    REQUIRE(cols % block_size == 0);
+    int const n_blocks = rows/block_size;
+    auto const A_block = vtint::split_to_block_matrix(A_h, n_blocks,n_blocks);
+    auto const b_block = vtint::split_to_block_vector(b_h, n_blocks);
+
+    //get the solution via solving A_h u = b_h
+    Eigen::SparseLU<Eigen::SparseMatrix<double>> solver;
+    solver.compute(A_h);
+    Eigen::VectorXd const sol_direct = solver.solve(b_h);
+
+
+    SECTION("Gauss Seidel"){
+        auto u0 = vtint::remove_initial_condition(u_initial);
+        u0 *= 0.;
+        for(int i = 0; i < 1; ++i){
+            u0 = vtint::gauss_seidel_iteration(A_block, b_block, u0);
+        }
+        SECTION("After one iteration, the equation system should be solved"){
+            auto const sol_smoother = u0.to_vector();
+            auto const diff = sol_direct - sol_smoother;
+            double const diff_norm = diff.lpNorm<Eigen::Infinity>();
+            REQUIRE(diff_norm < 1e-10);
+        }
+        REQUIRE(1==1);
+    }
+    SECTION("Jacobi"){
+        auto u0 = vtint::remove_initial_condition(u_initial);
+        u0 *= 0.;
+        SECTION("After one iteration, the equation system should not solved"){
+
+            for(int i = 0; i < 1; ++i){
+                u0 = vtint::jacobi_iteration(A_block, b_block, u0);
+            }
+            auto const sol_smoother = u0.to_vector();
+            auto const diff = sol_direct - sol_smoother;
+            double const diff_norm = diff.lpNorm<Eigen::Infinity>();
+            REQUIRE(diff_norm > 1e-3);
+        }
+        SECTION("After some iterations, the equation system should be approximately solved"){
+
+            for(int i = 0; i < 10; ++i){
+                u0 = vtint::jacobi_iteration(A_block, b_block, u0);
+            }
+            auto const sol_smoother = u0.to_vector();
+            auto const diff = sol_direct - sol_smoother;
+            double const diff_norm = diff.lpNorm<Eigen::Infinity>();
+            REQUIRE(diff_norm < 1e-10);
+        }
+        SECTION("After some iterations, the equation system should be approximately solved by relaxed scheme"){
+            double constexpr omega = .8;
+            for(int i = 0; i < 20; ++i){
+                u0 = vtint::relaxed_smoother_iteration(omega, vtint::jacobi_iteration, A_block, b_block, u0);
+            }
+            auto const sol_smoother = u0.to_vector();
+            auto const diff = sol_direct - sol_smoother;
+            double const diff_norm = diff.lpNorm<Eigen::Infinity>();
+            REQUIRE(diff_norm < 5e-10);
+        }
+        REQUIRE(1==1);
+    }
+
+
+}
\ No newline at end of file

simple_multigrid.cpp

+//
+// Created by mreic on 04.09.2024.
+//
+
+#include "simple_multigrid.hpp"
+#include <Eigen/SparseLU>
+
+namespace vtint {
+    grid_hierarchy_t::grid_hierarchy_t(const int n_levels, mfem::FiniteElementSpace &T_ansatz,
+                                       mfem::FiniteElementSpace &T_test) {
+        if (T_ansatz.GetMesh() != T_test.GetMesh()) {
+            throw std::runtime_error("T_ansatz and T_test must have the same mesh");
+        }
+        levels.resize(n_levels);
+        mfem::Mesh *mesh_ptr = T_ansatz.GetMesh();
+        T_ansatz_for_level.push_back(T_ansatz);
+        T_test_for_level.push_back(T_test);
+        for (int i = 1; i < n_levels; ++i) {
+            mesh_ptr->UniformRefinement();
+            auto T_ansatz_new = mfem::FiniteElementSpace(mesh_ptr, T_ansatz.FEColl());
+            auto T_test_new = mfem::FiniteElementSpace(mesh_ptr, T_test.FEColl());
+            T_ansatz_for_level.push_back(T_ansatz_new);
+            T_test_for_level.push_back(T_test_new);
+
+            mfem::OperatorHandle P_ansatz;
+            T_ansatz_new.GetTransferOperator(T_ansatz_for_level[i-1],P_ansatz);
+
+            mfem::OperatorHandle P_test;
+            T_test_new.GetTransferOperator(T_test_for_level[i-1],P_test);
+
+            levels[i-1].P = vtint::mfem_to_eigen(*P_ansatz.As<mfem::SparseMatrix>());
+            levels[i].R = vtint::mfem_to_eigen(*P_ansatz.As<mfem::SparseMatrix>()).transpose();
+
+        }
+
+        int tmpdbg = 0;
+    }
+
+
+    vtint::block_vector_t gauss_seidel_iteration(vtint::block_matrix_t const& A, vtint::block_vector_t const& b, vtint::block_vector_t const& u0){
+        vtint::block_vector_t u = u0;
+        for(int i=0; i<A.get_n_rows(); ++i){
+            Eigen::SparseMatrix<double> const &Aii = A.get_block(i,i);
+            Eigen::VectorXd const &bi = b.get_block(i);
+            Eigen::VectorXd bi_new = bi;
+            for(int j=0; j<A.get_n_cols(); ++j){
+                if(j!=i){
+                    bi_new -= A.get_block(i,j)*u.get_block(j);
+                }
+            }
+            //solve spatial system with sparse LU
+            Eigen::SparseLU<Eigen::SparseMatrix<double>> solver;
+            solver.compute(Aii);
+            u.get_block(i) = solver.solve(bi_new);
+        }
+        return u;
+    }
+
+
+    vtint::block_vector_t jacobi_iteration(vtint::block_matrix_t const& A, vtint::block_vector_t const& b, vtint::block_vector_t const& u0){
+        vtint::block_vector_t u = u0;
+        for(int i=0; i<A.get_n_rows(); ++i){
+            Eigen::SparseMatrix<double> const &Aii = A.get_block(i,i);
+            Eigen::VectorXd const &bi = b.get_block(i);
+            Eigen::VectorXd bi_new = bi;
+            for(int j=0; j<A.get_n_cols(); ++j){
+                if(j!=i){
+                    bi_new -= A.get_block(i,j)*u0.get_block(j);
+                }
+            }
+            //solve spatial system with sparse LU
+            Eigen::SparseLU<Eigen::SparseMatrix<double>> solver;
+            solver.compute(Aii);
+            u.get_block(i) = solver.solve(bi_new);
+        }
+        return u;
+    }
+
+
+
+} // vtint
\ No newline at end of file

simple_multigrid.hpp

+//
+// Created by mreic on 04.09.2024.
+//
+
+#ifndef VTINT_SIMPLE_MULTIGRID_HPP
+#define VTINT_SIMPLE_MULTIGRID_HPP
+
+#include "block_matrix_t.hpp"
+#include "block_vector_t.hpp"
+#include "eigen_helper.hpp"
+#include "mfem.hpp"
+
+#include <memory>
+
+
+namespace vtint {
+
+    struct level_t {
+        block_matrix_t A;
+        Eigen::SparseMatrix<double> P; //prolongation to finer level
+        Eigen::SparseMatrix<double> R; //restriction to coarser level
+    };
+
+
+    class grid_hierarchy_t{
+    public:
+        typedef std::shared_ptr<mfem::FiniteElementSpace> space_ptr_t;
+        grid_hierarchy_t(int const n_levels, mfem::FiniteElementSpace& T_ansatz, mfem::FiniteElementSpace& T_test);
+    private:
+        std::vector<level_t> levels;
+        std::vector<mfem::FiniteElementSpace> T_ansatz_for_level;
+        std::vector<mfem::FiniteElementSpace> T_test_for_level;
+    };
+
+    vtint::block_vector_t gauss_seidel_iteration(vtint::block_matrix_t const& A, vtint::block_vector_t const& b, vtint::block_vector_t const& u0);
+    vtint::block_vector_t jacobi_iteration(vtint::block_matrix_t const& A, vtint::block_vector_t const& b, vtint::block_vector_t const& u0);
+
+    template <class S_t>
+    vtint::block_vector_t relaxed_smoother_iteration(double const omega, S_t const& smoother, vtint::block_matrix_t const& A, vtint::block_vector_t const& b, vtint::block_vector_t const& u0){
+        vtint::block_vector_t u1 = smoother(A,b,u0);
+        return (1-omega)*u0 + omega*u1;
+    }
+
+} // vtint
+
+#endif //VTINT_SIMPLE_MULTIGRID_HPP

vtint/CMakeLists.txt

@@ -6,7 +6,9 @@ add_library(vtint vtint.cpp vtint.hpp common.hpp domain_integrators.cpp domain_i
         time_dependent_gridfunction.cpp
         time_dependent_gridfunction.hpp
         time_dependent_transformed_gridcoefficient_matrix.cpp
-        time_dependent_transformed_gridcoefficient_matrix.hpp)
+        time_dependent_transformed_gridcoefficient_matrix.hpp
+        ../simple_multigrid.cpp
+        ../simple_multigrid.hpp)
 TARGET_LINK_LIBRARIES(vtint ${MFEM_LIBRARIES} block_structures eigen_helper mfem_helper nonlinear_solver)
 target_include_directories(vtint
         PUBLIC ${MFEM_INCLUDE_DIRS}