From 0a2fb2d44785d7c9fab31cbbf6012f0da3e6472a Mon Sep 17 00:00:00 2001 From: Saveliy Date: Sun, 26 May 2024 00:49:52 +0300 Subject: [PATCH 1/8] update accuracy --- Poisson_Equation/headers/Solvers.hpp | 28 ++++++++++++++++++++++++---- Poisson_Equation/src/test.cpp | 13 +++++++++---- 2 files changed, 33 insertions(+), 8 deletions(-) diff --git a/Poisson_Equation/headers/Solvers.hpp b/Poisson_Equation/headers/Solvers.hpp index 161c33b..a3ce002 100644 --- a/Poisson_Equation/headers/Solvers.hpp +++ b/Poisson_Equation/headers/Solvers.hpp @@ -60,6 +60,18 @@ namespace numcpp return norm; } + FP error(const std::vector& v1, const std::vector& 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 vector_FMA(const std::vector& lhs, FP coef, const std::vector& rhs) { const size_t size = lhs.size(); @@ -210,15 +222,19 @@ namespace numcpp for (size_t i = 0; i < max_iterations; ++i) { + std::vector saved_approximation = approximation; + std::vector residual = (*system_matrix) * approximation - b; - FP residual_norm = norm(residual); - if (residual_norm <= required_precision) break; + //FP residual_norm = norm(residual); std::vector 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; @@ -439,10 +455,11 @@ namespace numcpp for (size_t i = 1; i < max_iterations; ++i) { + std::vector saved_approximation = approximation; + residual = (*system_matrix) * approximation - b; - FP residual_norm = norm(residual); - if (residual_norm <= required_precision) break; + //FP residual_norm = norm(residual); FP beta = scalar_product(Ah, residual) / scalar_product(Ah, h); @@ -453,6 +470,9 @@ 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) break; } return approximation; diff --git a/Poisson_Equation/src/test.cpp b/Poisson_Equation/src/test.cpp index b20b95e..e3b72ea 100644 --- a/Poisson_Equation/src/test.cpp +++ b/Poisson_Equation/src/test.cpp @@ -120,10 +120,15 @@ void test_task_custom_grid(size_t m, size_t n, std::unique_ptr int main() { - size_t m = 1024; - size_t n = 1024; + size_t m1 = 1024; + size_t n1 = 1024; - auto LS_solver = std::make_unique(std::vector(), 1000000, 0.00001, nullptr, std::vector()); + auto LS_solver1 = std::make_unique(std::vector(), 1000000, 0.000000001, nullptr, std::vector()); + test_task_custom_grid(m1, n1, std::move(LS_solver1)); - test_task_custom_grid(m, n, std::move(LS_solver)); + size_t m2 = 200; + size_t n2 = 200; + + auto LS_solver2 = std::make_unique(std::vector(), 1000000, 0.000001, nullptr, std::vector()); + test_task(m2, n2, std::move(LS_solver2)); } From 9635d721f758b089e0c95d7ec022b68dd8386c69 Mon Sep 17 00:00:00 2001 From: BlackAngel2108 Date: Sun, 26 May 2024 18:07:51 +0300 Subject: [PATCH 2/8] chebyshev test --- Poisson_Equation/headers/Solvers.hpp | 34 +++++++++++------- Poisson_Equation/src/test.cpp | 52 +++++++++++++++++++++++----- 2 files changed, 66 insertions(+), 20 deletions(-) diff --git a/Poisson_Equation/headers/Solvers.hpp b/Poisson_Equation/headers/Solvers.hpp index a3ce002..dae1fa4 100644 --- a/Poisson_Equation/headers/Solvers.hpp +++ b/Poisson_Equation/headers/Solvers.hpp @@ -388,8 +388,8 @@ namespace numcpp FP Mmin; FP Mmax; public: - ChebyshevIteration(std::vector initial_approximation, size_t max_iterations, FP required_precision, std::unique_ptr system_matrix, std::vector b) : - ISolver(std::move(initial_approximation), max_iterations, required_precision, std::move(system_matrix), std::move(b)) {} + ChebyshevIteration(std::vector initial_approximation, size_t max_iterations, FP required_precision, std::unique_ptr system_matrix, std::vector 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; @@ -402,28 +402,38 @@ namespace numcpp std::vector solve() const override { std::vector 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 = "< residual = (*system_matrix) * approximation - b; - - FP residual_norm = norm(residual); - if (residual_norm <= required_precision) break; + std::vector 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 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; } }; diff --git a/Poisson_Equation/src/test.cpp b/Poisson_Equation/src/test.cpp index e3b72ea..c311c15 100644 --- a/Poisson_Equation/src/test.cpp +++ b/Poisson_Equation/src/test.cpp @@ -118,17 +118,53 @@ void test_task_custom_grid(size_t m, size_t n, std::unique_ptr 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 corners = { -1.0, -1.0, 1.0, 1.0 }; + std::vector 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 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(init_app, 1000000000, 0.0000000000001, nullptr, std::vector(), Mmin, Mmax)); + + auto [duration, solution] = estimate_time(dirichlet_task); +} + int main() { - size_t m1 = 1024; - size_t n1 = 1024; + // size_t m1 = 1024; + // size_t n1 = 1024; + + // auto LS_solver1 = std::make_unique(std::vector(), 1000000, 0.000000001, nullptr, std::vector()); + // test_task_custom_grid(m1, n1, std::move(LS_solver1)); - auto LS_solver1 = std::make_unique(std::vector(), 1000000, 0.000000001, nullptr, std::vector()); - test_task_custom_grid(m1, n1, std::move(LS_solver1)); + // size_t m2 = 200; + // size_t n2 = 200; - size_t m2 = 200; - size_t n2 = 200; + // auto LS_solver2 = std::make_unique(std::vector(), 1000000, 0.000001, nullptr, std::vector()); + // test_task(m2, n2, std::move(LS_solver2)); - auto LS_solver2 = std::make_unique(std::vector(), 1000000, 0.000001, nullptr, std::vector()); - test_task(m2, n2, std::move(LS_solver2)); + test_Chebyshev(1024,1024); } From 46202ca52e56e5855b2079b66d1a864f0e9214aa Mon Sep 17 00:00:00 2001 From: Saveliy Date: Sun, 26 May 2024 20:06:34 +0300 Subject: [PATCH 3/8] update solution --- .../headers/Dirichlet_Problem.hpp | 46 ++++++++++++++++++- 1 file changed, 45 insertions(+), 1 deletion(-) diff --git a/Poisson_Equation/headers/Dirichlet_Problem.hpp b/Poisson_Equation/headers/Dirichlet_Problem.hpp index ef8be4b..289c6fe 100644 --- a/Poisson_Equation/headers/Dirichlet_Problem.hpp +++ b/Poisson_Equation/headers/Dirichlet_Problem.hpp @@ -508,8 +508,52 @@ std::vector> const DirichletProblemSolver::so } std::cout << "General error: " << approximation_error << std::endl; + std::vector> res(n + 1, std::vector(m + 1, 0)); + + for (size_t j = m / 2; j <= m; ++j) + { + res[0][j] = mu1(start_y + static_cast(j) * k); + } + for (size_t j = 0; j <= m / 2; ++j) + { + res[n / 2][j] = mu2(start_y + static_cast(j) * k); + } + for (size_t j = 0; j <= m; ++j) + { + res[n][j] = mu3(start_y + static_cast(j) * k); + } + + for (size_t i = n / 2 + 1; i < n; ++i) + { + res[i][0] = mu4(start_x + static_cast(i) * h); + } + for (size_t i = 1; i < n / 2; ++i) + { + res[i][m / 2] = mu5(start_x + static_cast(i) * h); + } + for (size_t i = 1; i < n; ++i) + { + res[i][m] = mu6(start_x + static_cast(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>(); + return res; } From 3fdb3f33136ef4ab11b46f564c0a0db1e105b83a Mon Sep 17 00:00:00 2001 From: Saveliy Date: Sun, 26 May 2024 20:42:55 +0300 Subject: [PATCH 4/8] update test --- Poisson_Equation/src/test.cpp | 20 ++++++++++---------- 1 file changed, 10 insertions(+), 10 deletions(-) diff --git a/Poisson_Equation/src/test.cpp b/Poisson_Equation/src/test.cpp index c311c15..7e8b68b 100644 --- a/Poisson_Equation/src/test.cpp +++ b/Poisson_Equation/src/test.cpp @@ -154,17 +154,17 @@ void test_Chebyshev(size_t n, size_t m) { int main() { - // size_t m1 = 1024; - // size_t n1 = 1024; + size_t m = 200; + size_t n = 200; - // auto LS_solver1 = std::make_unique(std::vector(), 1000000, 0.000000001, nullptr, std::vector()); - // test_task_custom_grid(m1, n1, std::move(LS_solver1)); + std::cout << "Conjugate gradient method" << std::endl; + auto LS_solver1 = std::make_unique(std::vector(), 1000000, 0.000000001, nullptr, std::vector()); + test_task_custom_grid(m, n, std::move(LS_solver1)); - // size_t m2 = 200; - // size_t n2 = 200; + std::cout << "\nMinimal residual method" << std::endl; + auto LS_solver2 = std::make_unique(std::vector(), 1000000, 0.000001, nullptr, std::vector()); + test_task(m, n, std::move(LS_solver2)); - // auto LS_solver2 = std::make_unique(std::vector(), 1000000, 0.000001, nullptr, std::vector()); - // test_task(m2, n2, std::move(LS_solver2)); - - test_Chebyshev(1024,1024); + std::cout << "\nChebyshev iteration method" << std::endl; + test_Chebyshev(m, n); } From f839c0eeb24398fe6e3996593b672a206f9d0495 Mon Sep 17 00:00:00 2001 From: Saveliy Date: Mon, 27 May 2024 15:11:24 +0300 Subject: [PATCH 5/8] test mcg --- Poisson_Equation/headers/Dirichlet_Problem.hpp | 18 +++++++++++++++++- Poisson_Equation/headers/Solvers.hpp | 10 ++++++++-- Poisson_Equation/src/test.cpp | 16 ++++++++-------- 3 files changed, 33 insertions(+), 11 deletions(-) diff --git a/Poisson_Equation/headers/Dirichlet_Problem.hpp b/Poisson_Equation/headers/Dirichlet_Problem.hpp index 289c6fe..3878819 100644 --- a/Poisson_Equation/headers/Dirichlet_Problem.hpp +++ b/Poisson_Equation/headers/Dirichlet_Problem.hpp @@ -501,11 +501,27 @@ std::vector> const DirichletProblemSolver::so solver->set_b(b); std::vector 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(n / 2 + 1 + i % (n / 2 - 1)) * h; + y_max_err = start_y + static_cast(1 + static_cast(i / (n / 2 - 1))) * k; + } + else + { + x_max_err = start_x + static_cast(1 + (i - (n / 2 - 1) * (m / 2)) % (n - 1)) * h; + y_max_err = start_y + static_cast(m / 2 + 1 + static_cast((i - (n / 2 - 1)) / (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> res(n + 1, std::vector(m + 1, 0)); diff --git a/Poisson_Equation/headers/Solvers.hpp b/Poisson_Equation/headers/Solvers.hpp index dae1fa4..9ea3fe4 100644 --- a/Poisson_Equation/headers/Solvers.hpp +++ b/Poisson_Equation/headers/Solvers.hpp @@ -465,11 +465,12 @@ namespace numcpp for (size_t i = 1; i < max_iterations; ++i) { + std::cout << "======= ITER " << i << " =======" << std::endl; std::vector saved_approximation = approximation; residual = (*system_matrix) * approximation - b; - //FP residual_norm = norm(residual); + FP residual_norm = norm(residual); FP beta = scalar_product(Ah, residual) / scalar_product(Ah, h); @@ -482,7 +483,12 @@ namespace numcpp approximation = vector_FMA(h, alpha, approximation); FP approximation_error = error(saved_approximation, approximation); - if (approximation_error <= required_precision) break; + if (approximation_error <= required_precision) + { + std::cout << "Approximation error: " << approximation_error << std::endl; + std::cout << "Residual norm: " << residual_norm << std::endl; + break; + } } return approximation; diff --git a/Poisson_Equation/src/test.cpp b/Poisson_Equation/src/test.cpp index 7e8b68b..27b7676 100644 --- a/Poisson_Equation/src/test.cpp +++ b/Poisson_Equation/src/test.cpp @@ -154,17 +154,17 @@ void test_Chebyshev(size_t n, size_t m) { int main() { - size_t m = 200; - size_t n = 200; + size_t m = 2048; + size_t n = 2048; std::cout << "Conjugate gradient method" << std::endl; - auto LS_solver1 = std::make_unique(std::vector(), 1000000, 0.000000001, nullptr, std::vector()); + auto LS_solver1 = std::make_unique(std::vector(), 1000000, 0.0000000000001, nullptr, std::vector()); test_task_custom_grid(m, n, std::move(LS_solver1)); - std::cout << "\nMinimal residual method" << std::endl; - auto LS_solver2 = std::make_unique(std::vector(), 1000000, 0.000001, nullptr, std::vector()); - test_task(m, n, std::move(LS_solver2)); + //std::cout << "\nMinimal residual method" << std::endl; + //auto LS_solver2 = std::make_unique(std::vector(), 1000000, 0.000001, nullptr, std::vector()); + //test_task(m, n, std::move(LS_solver2)); - std::cout << "\nChebyshev iteration method" << std::endl; - test_Chebyshev(m, n); + //std::cout << "\nChebyshev iteration method" << std::endl; + //test_Chebyshev(m, n); } From c566d83f913e4faa1647025c2d6edfbe56a161b4 Mon Sep 17 00:00:00 2001 From: Saveliy Borisov <94908850+Amazingkivas@users.noreply.github.com> Date: Mon, 27 May 2024 15:53:53 +0300 Subject: [PATCH 6/8] Update test.cpp --- Poisson_Equation/src/test.cpp | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/Poisson_Equation/src/test.cpp b/Poisson_Equation/src/test.cpp index 27b7676..744640e 100644 --- a/Poisson_Equation/src/test.cpp +++ b/Poisson_Equation/src/test.cpp @@ -154,8 +154,8 @@ void test_Chebyshev(size_t n, size_t m) { int main() { - size_t m = 2048; - size_t n = 2048; + size_t m = 4096; + size_t n = 4096; std::cout << "Conjugate gradient method" << std::endl; auto LS_solver1 = std::make_unique(std::vector(), 1000000, 0.0000000000001, nullptr, std::vector()); From 9b663972b02e5b7f1b80f9feee695634e0988bec Mon Sep 17 00:00:00 2001 From: Saveliy Borisov <94908850+Amazingkivas@users.noreply.github.com> Date: Wed, 29 May 2024 00:25:59 +0300 Subject: [PATCH 7/8] Update test.cpp --- Poisson_Equation/src/test.cpp | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/Poisson_Equation/src/test.cpp b/Poisson_Equation/src/test.cpp index 744640e..27b7676 100644 --- a/Poisson_Equation/src/test.cpp +++ b/Poisson_Equation/src/test.cpp @@ -154,8 +154,8 @@ void test_Chebyshev(size_t n, size_t m) { int main() { - size_t m = 4096; - size_t n = 4096; + size_t m = 2048; + size_t n = 2048; std::cout << "Conjugate gradient method" << std::endl; auto LS_solver1 = std::make_unique(std::vector(), 1000000, 0.0000000000001, nullptr, std::vector()); From fb7e316f224dfa996f5efc7dc1fccf7e4a4c555f Mon Sep 17 00:00:00 2001 From: Saveliy Borisov <94908850+Amazingkivas@users.noreply.github.com> Date: Wed, 29 May 2024 00:27:35 +0300 Subject: [PATCH 8/8] Update Dirichlet_Problem.hpp --- Poisson_Equation/headers/Dirichlet_Problem.hpp | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/Poisson_Equation/headers/Dirichlet_Problem.hpp b/Poisson_Equation/headers/Dirichlet_Problem.hpp index 3878819..0b0973f 100644 --- a/Poisson_Equation/headers/Dirichlet_Problem.hpp +++ b/Poisson_Equation/headers/Dirichlet_Problem.hpp @@ -517,7 +517,7 @@ std::vector> const DirichletProblemSolver::so else { x_max_err = start_x + static_cast(1 + (i - (n / 2 - 1) * (m / 2)) % (n - 1)) * h; - y_max_err = start_y + static_cast(m / 2 + 1 + static_cast((i - (n / 2 - 1)) / (n - 1))) * k; + y_max_err = start_y + static_cast(m / 2 + 1 + static_cast((i - (n / 2 - 1) * (m / 2)) / (n - 1))) * k; } } }