Test mcg

Poisson_Equation/headers/Dirichlet_Problem.hpp

@@ -501,15 +501,75 @@ std::vector<std::vector<FP>> const DirichletProblemSolver<GridType::Regular>::so
         solver->set_b(b);
         std::vector<FP> solution = solver->solve();
 
+        FP x_max_err = 0.0;
+        FP y_max_err = 0.0;
         FP approximation_error = 0.0;
         for (size_t i = 0; i < sz; ++i)
         {
-            approximation_error = std::max(std::abs(solution[i] - real_sol[i]), approximation_error);
+            if (std::abs(solution[i] - real_sol[i]) > approximation_error)
+            {
+                approximation_error = std::abs(solution[i] - real_sol[i]);
+                if (i < (n / 2 - 1) * (m / 2))
+                {
+                    x_max_err = start_x + static_cast<FP>(n / 2 + 1 + i % (n / 2 - 1)) * h;
+                    y_max_err = start_y + static_cast<FP>(1 + static_cast<int>(i / (n / 2 - 1))) * k;
+                }
+                else
+                {
+                    x_max_err = start_x + static_cast<FP>(1 + (i - (n / 2 - 1) * (m / 2)) % (n - 1)) * h;
+                    y_max_err = start_y + static_cast<FP>(m / 2 + 1 + static_cast<int>((i - (n / 2 - 1) * (m / 2)) / (n - 1))) * k;
+                }
+            }
         }
+        std::cout << "Node with maximum error: x = " << x_max_err << ", y = " << y_max_err << std::endl;
         std::cout << "General error: " << approximation_error << std::endl;
 
+        std::vector<std::vector<FP>> res(n + 1, std::vector<FP>(m + 1, 0));
+
+        for (size_t j = m / 2; j <= m; ++j)
+        {
+            res[0][j] = mu1(start_y + static_cast<FP>(j) * k);
+        }
+        for (size_t j = 0; j <= m / 2; ++j)
+        {
+            res[n / 2][j] = mu2(start_y + static_cast<FP>(j) * k);
+        }
+        for (size_t j = 0; j <= m; ++j)
+        {
+            res[n][j] = mu3(start_y + static_cast<FP>(j) * k);
+        }
+
+        for (size_t i = n / 2 + 1; i < n; ++i)
+        {
+            res[i][0] = mu4(start_x + static_cast<FP>(i) * h);
+        }
+        for (size_t i = 1; i < n / 2; ++i)
+        {
+            res[i][m / 2] = mu5(start_x + static_cast<FP>(i) * h);
+        }
+        for (size_t i = 1; i < n; ++i)
+        {
+            res[i][m] = mu6(start_x + static_cast<FP>(i) * h);
+        }
+
+        size_t index = 0;
+        for (size_t j = 1; j <= m / 2; ++j)
+        {
+            for (size_t i = n / 2 + 1; i < n; ++i)
+            {
+                res[i][j] = solution[index++];
+            }
+        }
+        for (size_t j = m / 2 + 1; j < m; ++j)
+        {
+            for (size_t i = 1; i < n; ++i)
+            {
+                res[i][j] = solution[index++];
+            }
+        }
+
         // Placeholder
-        return std::vector<std::vector<FP>>();
+        return res;
     }
 
 

Poisson_Equation/headers/Solvers.hpp

@@ -60,6 +60,18 @@ namespace numcpp
         return norm;
     }
 
+    FP error(const std::vector<FP>& v1, const std::vector<FP>& v2)
+    {
+        FP error = 0.0;
+
+#pragma omp parallel for reduction(max: error)
+        for (size_t i = 0; i < v1.size(); ++i)
+        {
+            error = std::max(std::abs(v1[i] - v2[i]), error);
+        }
+        return error;
+    }
+
     std::vector<FP> vector_FMA(const std::vector<FP>& lhs, FP coef, const std::vector<FP>& rhs)
     {
         const size_t size = lhs.size();
@@ -210,15 +222,19 @@ namespace numcpp
 
             for (size_t i = 0; i < max_iterations; ++i)
             {
+                std::vector<FP> saved_approximation = approximation;
+
                 std::vector<FP> residual = (*system_matrix) * approximation - b;
 
-                FP residual_norm = norm(residual);
-                if (residual_norm <= required_precision) break;
+                //FP residual_norm = norm(residual);
 
                 std::vector<FP> Ar = (*system_matrix) * residual;
                 FP tau = scalar_product(Ar, residual) / scalar_product(Ar, Ar);
 
                 approximation = vector_FMA(residual, -tau, approximation);
+
+                FP approximation_error = error(saved_approximation, approximation);
+                if (approximation_error <= required_precision) break;
             }
 
             return approximation;
@@ -372,8 +388,8 @@ namespace numcpp
         FP Mmin;
         FP Mmax;
     public:
-        ChebyshevIteration(std::vector<FP> initial_approximation, size_t max_iterations, FP required_precision, std::unique_ptr<IMatrix> system_matrix, std::vector<FP> b) :
-            ISolver(std::move(initial_approximation), max_iterations, required_precision, std::move(system_matrix), std::move(b)) {}
+        ChebyshevIteration(std::vector<FP> initial_approximation, size_t max_iterations, FP required_precision, std::unique_ptr<IMatrix> system_matrix, std::vector<FP> b, FP min = 0, FP max = 0) :
+            ISolver(std::move(initial_approximation), max_iterations, required_precision, std::move(system_matrix), std::move(b)), Mmin(min), Mmax(max) {}
 
         ChebyshevIteration() = default;
         ~ChebyshevIteration() = default;
@@ -386,28 +402,38 @@ namespace numcpp
         std::vector<FP> solve() const override
         {
             std::vector<FP> approximation = std::move(initial_approximation);
-            double size = (*system_matrix).size();
-            //FP h = sqrt(1.0 / (*system_matrix).at(0, 1));
-            //FP k = sqrt(1.0 / (*system_matrix).at(0, sqrt(size)));
 
             FP k_cheb = 2.0;
             FP tau0 = 1.0 / ((Mmin + Mmax) / 2.0 + (Mmax - Mmin) / 2 * cos(PI / (2.0 * k_cheb) * (1.0 + 2.0 * 0.0)));
             FP tau1 = 1.0 / ((Mmin + Mmax) / 2.0 + (Mmax - Mmin) / 2 * cos(PI / (2.0 * k_cheb) * (1.0 + 2.0 * 1.0)));
 
-
+            std::cout<<"tau1 = "<<tau0<<" tau2 = "<<tau1<<"\n";
             for (size_t i = 0; i < max_iterations; ++i)
             {
                 std::vector<FP> residual = (*system_matrix) * approximation - b;
-
-                FP residual_norm = norm(residual);
-                if (residual_norm <= required_precision) break;
+                std::vector<FP> saved_approximation = approximation;
 
                 if (i % 2 == 0)
-                    approximation = vector_FMA(residual, -tau0, approximation);
+                    approximation = vector_FMA(residual, tau0, approximation);
                 else
-                    approximation = vector_FMA(residual, -tau1, approximation);
+                    approximation = vector_FMA(residual, tau1, approximation);
+
+                FP approximation_error = error(saved_approximation, approximation);
+                if (approximation_error <= required_precision) break;
             }
-
+            // for (size_t i = 0; i < max_iterations; ++i)
+            // {
+            //     std::vector<FP> residual = (*system_matrix) * approximation - b;
+
+            //     FP residual_norm = norm(residual);
+            //     if (residual_norm <= required_precision) break;
+
+            //     if (i % 2 == 0)
+            //         approximation = vector_FMA(residual, -tau0, approximation);
+            //     else
+            //         approximation = vector_FMA(residual, -tau1, approximation);
+            //     std::cout << residual_norm<< "  "<<  required_precision <<"\n";
+            // }
             return approximation;
         }
     };
@@ -439,10 +465,12 @@ namespace numcpp
 
             for (size_t i = 1; i < max_iterations; ++i)
             {
+                std::cout << "======= ITER " << i << " =======" << std::endl;
+                std::vector<FP> saved_approximation = approximation;
+
                 residual = (*system_matrix) * approximation - b;
 
                 FP residual_norm = norm(residual);
-                if (residual_norm <= required_precision) break;
 
                 FP beta = scalar_product(Ah, residual) / scalar_product(Ah, h);
 
@@ -453,6 +481,14 @@ namespace numcpp
                 alpha = -scalar_product(residual, h) / scalar_product(Ah, h);
 
                 approximation = vector_FMA(h, alpha, approximation);
+
+                FP approximation_error = error(saved_approximation, approximation);
+                if (approximation_error <= required_precision) 
+                {
+                    std::cout << "Approximation error: " << approximation_error << std::endl;
+                    std::cout << "Residual norm: " << residual_norm << std::endl;
+                    break;
+                }
             }
 
             return approximation;

Poisson_Equation/src/test.cpp

@@ -118,12 +118,53 @@ void test_task_custom_grid(size_t m, size_t n, std::unique_ptr<numcpp::ISolver>
     auto [duration, solution] = estimate_time(solver);
 }
 
+void test_Chebyshev(size_t n, size_t m) {
+    // size_t n = 2048;
+    // size_t m = 2048;
+    std::cout << "n = " << n << " m = " << m << '\n';
+
+    std::array<double, 4> corners = { -1.0, -1.0, 1.0, 1.0 };
+    std::vector<FP> init_app((n - 1) * (m - 1), 0.0);
+
+    FP h = (corners[2] - corners[0]) / n;
+    FP k = (corners[3] - corners[1]) / m;
+
+    FP Mmin = 4.0 / pow(h, 2) * pow(sin(numcpp::PI / 2.0 / n), 2) + 4.0 / pow(k, 2) * pow(sin(numcpp::PI / 2.0 / m), 2);
+    FP Mmax = 4.0 / pow(h, 2) * pow(sin(numcpp::PI * (n - 1) / 2.0 / n), 2) + 4.0 / pow(k, 2) * pow(sin(numcpp::PI * (m - 1) / 2.0 / m), 2);
+    std::cout << "Mmin = " << Mmin << " Mmax = " << Mmax << '\n';
+    auto u = [](double x, double y) { return exp(1 - pow(x, 2) - pow(y, 2)); };
+    auto f = [](double x, double y) { return -4 * exp(1 - pow(x, 2) - pow(y, 2)) * (pow(y, 2) + pow(x, 2) - 1); };
+
+    auto mu1 = [](double y) { return exp(-pow(y, 2)); };
+    auto mu2 = [](double y) { return exp(-pow(y, 2)); };
+    auto mu3 = [](double x) { return exp(-pow(x, 2)); };
+    auto mu4 = [](double x) { return exp(-pow(x, 2)); };
+
+    numcpp::DirichletProblemSolver<numcpp::Regular> dirichlet_task;
+    dirichlet_task.set_fraction(m, n);
+    dirichlet_task.set_corners(corners);
+    dirichlet_task.set_u(u);
+    dirichlet_task.set_f(f);
+    dirichlet_task.set_boundary_conditions({ mu1, mu2, mu3, mu4 });
+
+    dirichlet_task.set_solver(std::make_unique<numcpp::ChebyshevIteration>(init_app, 1000000000, 0.0000000000001, nullptr, std::vector<FP>(), Mmin, Mmax));
+
+    auto [duration, solution] = estimate_time(dirichlet_task);
+}
+
 int main() 
 {
-    size_t m = 1024;
-    size_t n = 1024;
+    size_t m = 2048;
+    size_t n = 2048;
+
+    std::cout << "Conjugate gradient method" << std::endl;
+    auto LS_solver1 = std::make_unique<numcpp::ConGrad>(std::vector<FP>(), 1000000, 0.0000000000001, nullptr, std::vector<FP>());
+    test_task_custom_grid(m, n, std::move(LS_solver1));
 
-    auto LS_solver = std::make_unique<numcpp::ConGrad>(std::vector<FP>(), 1000000, 0.00001, nullptr, std::vector<FP>()); 
+    //std::cout << "\nMinimal residual method" << std::endl;
+    //auto LS_solver2 = std::make_unique<numcpp::MinRes>(std::vector<FP>(), 1000000, 0.000001, nullptr, std::vector<FP>());
+    //test_task(m, n, std::move(LS_solver2));
 
-    test_task_custom_grid(m, n, std::move(LS_solver));
+    //std::cout << "\nChebyshev iteration method" << std::endl;
+    //test_Chebyshev(m, n);
 }