diff --git a/.github/workflows/main.yml b/.github/workflows/main.yml new file mode 100644 index 0000000..254bf35 --- /dev/null +++ b/.github/workflows/main.yml @@ -0,0 +1,20 @@ +name: Build application + +on: [push, pull_request] + +jobs: + ubuntu-gcc-build: + runs-on: ubuntu-latest + steps: + - uses: actions/checkout@v1 + - name: Build + run: | + sudo apt-get update + sudo apt install libomp-dev + cd Poisson_Equation/src + g++ -mavx2 -mfma -fopenmp -o out.exe test.cpp -I./../headers/ -std=c++17 -O3 + - name: Run test + shell: bash + run: | + cd Poisson_Equation/src + ./out.exe diff --git a/Poisson_Equation/files/Approximation.txt b/Poisson_Equation/files/Approximation.txt new file mode 100644 index 0000000..e5455eb --- /dev/null +++ b/Poisson_Equation/files/Approximation.txt @@ -0,0 +1,10 @@ +0.367879 0.527292 0.697676 0.852144 0.960789 1 0.960789 0.852144 0.697676 0.527292 0.367879 +0.527292 0.755035 0.999793 1.2224 1.37932 1.43613 1.37967 1.22293 1.00029 0.755327 0.527292 +0.697676 0.999622 1.32406 1.61925 1.82745 1.90297 1.82831 1.62053 1.32523 1.00029 0.697676 +0.852144 1.22171 1.61832 1.97897 2.23333 2.3258 2.23516 1.98153 1.62053 1.22293 0.852144 +0.960789 1.37769 1.82441 2.23022 2.51625 2.6204 2.52012 2.23516 1.82831 1.37967 0.960789 +1 1.43333 1.89648 2.31637 2.6117 2.71828 2.6204 2.3258 1.90297 1.43613 1 +0 0 0 0 0 2.6117 2.51625 2.23333 1.82745 1.37932 0.960789 +0 0 0 0 0 2.31637 2.23022 1.97897 1.61925 1.2224 0.852144 +0 0 0 0 0 1.89648 1.82441 1.61832 1.32406 0.999793 0.697676 +0 0 0 0 0 1.43333 1.37769 1.22171 0.999622 0.755035 0.527292 diff --git a/Poisson_Equation/files/Correct.txt b/Poisson_Equation/files/Correct.txt new file mode 100644 index 0000000..4dc95b4 --- /dev/null +++ b/Poisson_Equation/files/Correct.txt @@ -0,0 +1,10 @@ +0.367879 0.527292 0.697676 0.852144 0.960789 1 0.960789 0.852144 0.697676 0.527292 0.367879 +0.527292 0.755784 1 1.2214 1.37713 1.43333 1.37713 1.2214 1 0.755784 0.527292 +0.697676 1 1.32313 1.61607 1.82212 1.89648 1.82212 1.61607 1.32313 1 0.697676 +0.852144 1.2214 1.61607 1.97388 2.22554 2.31637 2.22554 1.97388 1.61607 1.2214 0.852144 +0.960789 1.37713 1.82212 2.22554 2.50929 2.6117 2.50929 2.22554 1.82212 1.37713 0.960789 +1 1.43333 1.89648 2.31637 2.6117 2.71828 2.6117 2.31637 1.89648 1.43333 1 +0 0 0 0 0 2.6117 2.50929 2.22554 1.82212 1.37713 0.960789 +0 0 0 0 0 2.31637 2.22554 1.97388 1.61607 1.2214 0.852144 +0 0 0 0 0 1.89648 1.82212 1.61607 1.32313 1 0.697676 +0 0 0 0 0 1.43333 1.37713 1.2214 1 0.755784 0.527292 diff --git a/Poisson_Equation/files/Error.txt b/Poisson_Equation/files/Error.txt new file mode 100644 index 0000000..7f7e4d1 --- /dev/null +++ b/Poisson_Equation/files/Error.txt @@ -0,0 +1,2 @@ +Задача решена с погрешностью 0.0108329 +Максимальное отклонение точного и численного решений в точке x = 0.2, y = 0.6 diff --git a/Poisson_Equation/files/Solver_results.txt b/Poisson_Equation/files/Solver_results.txt new file mode 100644 index 0000000..67e45cf --- /dev/null +++ b/Poisson_Equation/files/Solver_results.txt @@ -0,0 +1,3 @@ +2.8212e-09 +2.49716e-08 +22 diff --git a/Poisson_Equation/files/Solver_results_2N.txt b/Poisson_Equation/files/Solver_results_2N.txt new file mode 100644 index 0000000..47c82e4 --- /dev/null +++ b/Poisson_Equation/files/Solver_results_2N.txt @@ -0,0 +1,3 @@ +9.71618e-09 +5.01926e-06 +582 diff --git a/Poisson_Equation/headers/Dirichlet_Problem.hpp b/Poisson_Equation/headers/Dirichlet_Problem.hpp index ef8be4b..04d621c 100644 --- a/Poisson_Equation/headers/Dirichlet_Problem.hpp +++ b/Poisson_Equation/headers/Dirichlet_Problem.hpp @@ -7,6 +7,7 @@ #include #include #include +#include using FP = double; @@ -53,6 +54,7 @@ namespace numcpp { using two_dim_function = std::function; size_t m, n; + size_t task_num; std::unique_ptr solver; std::array corners; // {x0, y0, xn, ym} @@ -105,6 +107,11 @@ namespace numcpp { this->mu5 = std::move(new_mu[4]); this->mu6 = std::move(new_mu[5]); } + + void set_task_num(size_t num) + { + this->task_num = num; + } }; @@ -210,23 +217,69 @@ std::vector> const DirichletProblemSolver::so solver->set_system_matrix(std::move(matrix)); solver->set_b(b); std::vector solution = solver->solve(); + std::vector> res(n + 1, std::vector(m + 1, 0)); + if(task_num == 0) + { + FP error = 0; + FP x_max = 0; + FP y_max = 0; - std::vector> res; - res.resize(m - 1); + for(size_t i = 0; i < m - 1; i++) + { + for(size_t j = 0; j < n - 1; j++) + { + FP y = corners[1] + (i + 1) * k; + FP x = corners[0] + (j + 1) * h; + FP dot_solution = solution[(n - 1) * i + j]; + FP dot_correct = u(x, y); + if(std::abs(dot_correct - dot_solution) > error) + { + error = std::abs(dot_correct - dot_solution); + x_max = x; + y_max = y; + } + } + } + std::ofstream fout_res("../files/Error.txt"); + fout_res << "Задача решена с погрешностью " << error << "\nМаксимальное отклонение точного и численного решений в точке x = " << x_max << ", y = " << y_max << '\n'; + fout_res.close(); + } + + + + FP start_x = corners[0]; + FP start_y = corners[1]; - FP approximation_error = 0; + for (size_t j = 0; j <= m; ++j) + { + res[0][j] = mu1(start_y + static_cast(j) * k); + } + for (size_t j = 0; j <= m; ++j) + { + res[n][j] = mu2(start_y + static_cast(j) * k); + } - for(size_t i = 0; i < m - 1; i++) - for(size_t j = 0; j < n - 1; j++) + + for (size_t i = 1; i <= n; ++i) + { + res[i][0] = mu3(start_x + static_cast(i) * h); + } + for (size_t i = 1; i <= n; ++i) + { + res[i][m] = mu4(start_x + static_cast(i) * h); + } + + size_t index = 0; + for (size_t j = 1; j < m; ++j) + { + for (size_t i = 1; i < n; ++i) { - FP y = corners[1] + (i + 1) * k; - FP x = corners[0] + (j + 1) * h; - FP dot_solution = solution[(n - 1) * i + j]; - approximation_error = std::max(std::abs(u(x, y) - dot_solution), approximation_error); - res[i].push_back(dot_solution); + res[i][j] = solution[index++]; } - - std::cout << "General error: " << approximation_error << '\n'; + } + + + return res; } @@ -501,15 +554,80 @@ 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) * (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::ofstream fout_res("../files/Error.txt"); + fout_res << "Задача решена с погрешностью " << approximation_error << "\nМаксимальное отклонение точного и численного решений в точке x = " << x_max_err << ", y = " << y_max_err << '\n'; + fout_res.close(); + + 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; } @@ -517,4 +635,4 @@ std::vector> const DirichletProblemSolver::so -#endif // __DIRICHLET_PROBLEM_HPP__ +#endif // __DIRICHLET_PROBLEM_HPP__ \ No newline at end of file diff --git a/Poisson_Equation/headers/Solvers.hpp b/Poisson_Equation/headers/Solvers.hpp index 161c33b..a58d52a 100644 --- a/Poisson_Equation/headers/Solvers.hpp +++ b/Poisson_Equation/headers/Solvers.hpp @@ -8,6 +8,7 @@ #include #include #include +#include #include #include "./Dirichlet_Problem.hpp" @@ -60,6 +61,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(); @@ -156,6 +169,7 @@ namespace numcpp system_matrix(std::move(system_matrix)), b(b) {} + public: virtual std::vector solve() const = 0; @@ -197,6 +211,7 @@ namespace numcpp // Minimal residual method class MinRes : public ISolver { + public: MinRes(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)) {} @@ -207,20 +222,30 @@ namespace numcpp std::vector solve() const override { std::vector approximation = std::move(initial_approximation); + std::vector residual; + FP approximation_error, residual_norm; + size_t i = 0; - for (size_t i = 0; i < max_iterations; ++i) + for (; i < max_iterations; ++i) { - std::vector residual = (*system_matrix) * approximation - b; + std::vector saved_approximation = approximation; - FP residual_norm = norm(residual); - if (residual_norm <= required_precision) break; + residual = (*system_matrix) * approximation - b; std::vector Ar = (*system_matrix) * residual; FP tau = scalar_product(Ar, residual) / scalar_product(Ar, Ar); approximation = vector_FMA(residual, -tau, approximation); + + approximation_error = error(saved_approximation, approximation); + if (approximation_error <= required_precision) break; } + std::ofstream fout("../files/Solver_results.txt"); + fout << approximation_error << '\n' << norm(residual) << '\n' << i; + fout.close(); + + return approximation; } }; @@ -239,14 +264,14 @@ namespace numcpp { size_t size = initial_approximation.size(); std::vector approximation = std::move(initial_approximation); + size_t i = 0; + FP approximation_error = 0; + std::vector residual; //top relaxation main part - for (size_t i = 0; i < max_iterations; ++i) + for (; i < max_iterations; ++i) { - std::vector residual = (*system_matrix) * approximation - b; - - FP residual_norm = norm(residual); - if (residual_norm <= required_precision) break; + std::vector saved_approximation = approximation; for(size_t ii = 0; ii < size; ++ii) { @@ -262,10 +287,14 @@ namespace numcpp approx_new /= system_matrix->at(ii, ii); approximation[ii] = approx_new; } + approximation_error = error(saved_approximation, approximation); + if (approximation_error <= required_precision) break; } - std::vector residual = (*system_matrix) * approximation - b; - FP residual_norm = norm(residual); - std::cout << "погрешность решения СЛАУ " << residual_norm << std::endl; + residual = (*system_matrix) * approximation - b; + + std::ofstream fout("../files/Solver_results.txt"); + fout << approximation_error << '\n' << norm(residual) << '\n' << i; + fout.close(); return approximation; } @@ -284,6 +313,7 @@ namespace numcpp FP omega; size_t n, m; + public: TopRelaxationOptimizedForDirichletRegularGrid(std::vector initial_approximation, size_t max_iterations, FP required_precision, std::unique_ptr system_matrix, std::vector b, two_dim_function new_f, one_dim_function new_mu1, one_dim_function new_mu2, one_dim_function new_mu3, one_dim_function new_mu4, @@ -298,10 +328,8 @@ namespace numcpp { size_t size = initial_approximation.size(); std::vector approximation = std::move(initial_approximation); - - std::vector residual = (*system_matrix) * approximation - b; - FP residual_norm = norm(residual); - if (residual_norm <= required_precision) return approximation; + std::vector residual; + FP residual_norm; std::vector> approximation_matrix(n + 1, std::vector(m + 1)); FP h = (corners[2] - corners[0]) / n; @@ -347,19 +375,16 @@ namespace numcpp } if (residual_norm <= required_precision) break; } - if(ii == max_iterations) std::cout << "достигнуто максимальное количество шагов!\n"; - std::cout << "eps " << residual_norm << std::endl; - std::cout << "итераций " << ii << '\n'; - ind = 0; for(size_t i = 1; i < m; ++i) for(size_t j = 1; j < n; ++j) approximation[ind++] = approximation_matrix[j][i]; residual = (*system_matrix) * approximation - b; - residual_norm = norm(residual); - std::cout << "невязка " << residual_norm << '\n'; + std::ofstream fout("../files/Solver_results.txt"); + fout << residual_norm << '\n' << norm(residual) << '\n' << ii; + fout.close(); return approximation; } }; @@ -371,10 +396,10 @@ namespace numcpp private: 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)) {} + public: + 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; @@ -387,27 +412,39 @@ namespace numcpp { 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 approximation_error = 0; + size_t i = 0; 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 = "< 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); + // if (residual_norm <= required_precision) break; 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); + + approximation_error = error(saved_approximation, approximation); + if (approximation_error <= required_precision) break; } + std::vector residual = (*system_matrix) * approximation - b; + + std::ofstream fout("../files/Solver_results.txt"); + fout << approximation_error << '\n' << norm(residual) << '\n' << i; + fout.close(); + return approximation; } }; @@ -416,6 +453,7 @@ namespace numcpp // Conjugate gradient method class ConGrad : public ISolver { + public: ConGrad(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)) {} @@ -426,6 +464,7 @@ namespace numcpp std::vector solve() const override { std::vector approximation = std::move(initial_approximation); + FP approximation_error = 0; std::vector residual = (*system_matrix) * approximation - b; @@ -436,13 +475,15 @@ namespace numcpp FP alpha = -scalar_product(residual, h) / scalar_product(Ah, h); approximation = vector_FMA(h, alpha, approximation); + size_t i = 1; - for (size_t i = 1; i < max_iterations; ++i) + for (; 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,8 +494,17 @@ namespace numcpp alpha = -scalar_product(residual, h) / scalar_product(Ah, h); approximation = vector_FMA(h, alpha, approximation); + + approximation_error = error(saved_approximation, approximation); + if (approximation_error <= required_precision) break; } + residual = (*system_matrix) * approximation - b; + + std::ofstream fout("../files/Solver_results.txt"); + fout << approximation_error << '\n' << norm(residual) << '\n' << i; + fout.close(); + return approximation; } }; @@ -463,4 +513,4 @@ namespace numcpp } // namespace numcpp -#endif // __SOLVER_HPP__ +#endif // __SOLVER_HPP__ \ No newline at end of file diff --git a/Poisson_Equation/headers/tasks.hpp b/Poisson_Equation/headers/tasks.hpp new file mode 100644 index 0000000..6e97d3d --- /dev/null +++ b/Poisson_Equation/headers/tasks.hpp @@ -0,0 +1,288 @@ +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include + +#include "../headers/Dirichlet_Problem.hpp" +#include "../headers/Solvers.hpp" + +#define FP double + + +void test_task(size_t solver_num, size_t n, size_t m, size_t max_iterations, FP eps, FP omega){ + std::array corners = {-1.0, -1.0, 1.0, 1.0}; + std::vector init_app((n - 1) * (m - 1), 0.0); + + 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_task_num(0); + dirichlet_task.set_boundary_conditions({mu1, mu2, mu3, mu4}); + + 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); + + if(solver_num == 0) + dirichlet_task.set_solver(std::make_unique + (init_app, max_iterations, eps, nullptr, std::vector(), f, mu1, mu2, mu3, mu4, n, m, corners, omega)); + else if(solver_num == 1) + dirichlet_task.set_solver(std::make_unique(init_app, max_iterations, eps, nullptr, std::vector())); + else if (solver_num == 2) + dirichlet_task.set_solver(std::make_unique(init_app, max_iterations, eps, nullptr, std::vector(), Mmin, Mmax)); + else + dirichlet_task.set_solver(std::make_unique(init_app, max_iterations, eps, nullptr, std::vector())); + + auto start = std::chrono::high_resolution_clock::now(); + auto solution = dirichlet_task.solve(); + auto end = std::chrono::high_resolution_clock::now(); + + int64_t duration = std::chrono::duration_cast(end - start).count(); + std::cout << "Overall time: " << duration / 1000000 << "s\n"; + + std::ofstream fapp("../files/Approximation.txt"); + std::ofstream fcurr("../files/Correct.txt"); + for(size_t j = m; j > 0; j--) + { + for(size_t i = 0; i <= n; i++) + { + FP x = corners[0] + i * h; + FP y = corners[1] + j* k; + fapp << solution[i][j] <<' '; + fcurr << u(x, y) << ' '; + } + fapp << '\n'; + fcurr << '\n'; + } + fapp.close(); + fcurr.close(); + + +} + + +void test_custom_task(size_t solver_num, size_t n, size_t m, size_t max_iterations, FP eps, FP omega){ + std::array corners = {-1.0, -1.0, 1.0, 1.0}; + size_t size = (n / 2 - 1) * (m - 1) + (n / 2) * (m / 2 - 1); + std::vector init_app(size, 0.0); + + 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(1.0 - pow(y, 2)); }; + auto mu3 = [](double y) { return exp(-pow(y, 2)); }; + auto mu4 = [](double x) { return exp(-pow(x, 2)); }; + auto mu5 = [](double x) { return exp(1.0 - pow(x, 2)); }; + auto mu6 = [](double x) { return exp(-pow(x, 2)); }; + + numcpp::DirichletProblemSolver dirichlet_task; + std::array, 6> arr{mu1, mu2, mu3, mu4, mu5, mu6}; + + 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_for_r_shaped_grid(arr); + + 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); + + if(solver_num == 0) + dirichlet_task.set_solver(std::make_unique(init_app, max_iterations, eps, nullptr, std::vector(), omega)); + else if(solver_num == 1) + dirichlet_task.set_solver(std::make_unique(init_app, max_iterations, eps, nullptr, std::vector())); + else if (solver_num == 2) + dirichlet_task.set_solver(std::make_unique(init_app, max_iterations, eps, nullptr, std::vector(), Mmin, Mmax)); + else + dirichlet_task.set_solver(std::make_unique(init_app, max_iterations, eps, nullptr, std::vector())); + + auto start = std::chrono::high_resolution_clock::now(); + auto solution = dirichlet_task.solve(); + auto end = std::chrono::high_resolution_clock::now(); + + int64_t duration = std::chrono::duration_cast(end - start).count(); + std::cout << "Overall time: " << duration / 1000000 << "s\n"; + + std::ofstream fapp("../files/Approximation.txt"); + std::ofstream fcurr("../files/Correct.txt"); + + for(size_t j = m; j > 0; j--) + { + for(size_t i = 0; i <= n; i++) + { + FP x = corners[0] + i * h; + FP y = corners[1] + j * k; + fapp << solution[i][j] <<' '; + if(i < n / 2 && j < m / 2) + fcurr << 0.0 << ' '; + else + fcurr << u(x, y) << ' '; + } + fapp << '\n'; + fcurr << '\n'; + } + fapp.close(); + fcurr.close(); +} + + +void main_task(size_t solver_num, size_t n, size_t m, size_t max_iterations, FP eps, FP omega){ + std::array corners = {-1.0, -1.0, 1.0, 1.0}; + std::vector init_app((n - 1) * (m - 1), 0.0); + + auto f = [](double x, double y) { return std::abs(pow(sin(numcpp::PI*x*y), 3)); }; + + auto mu1 = [](double y) { return -1*pow(y, 2)+1; }; + auto mu2 = [](double y) { return -1*pow(y, 2)+1; }; + auto mu3 = [](double x) { return -1*pow(x, 2)+1; }; + auto mu4 = [](double x) { return -1*pow(x, 2)+1; }; + numcpp::DirichletProblemSolver dirichlet_task; + + dirichlet_task.set_fraction(m, n); + dirichlet_task.set_corners(corners); + dirichlet_task.set_f(f); + dirichlet_task.set_task_num(1); + dirichlet_task.set_boundary_conditions({mu1, mu2, mu3, mu4}); + + 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); + + if(solver_num == 0) + dirichlet_task.set_solver(std::make_unique + (init_app, max_iterations, eps, nullptr, std::vector(), f, mu1, mu2, mu3, mu4, n, m, corners, omega)); + else if(solver_num == 1) + dirichlet_task.set_solver(std::make_unique(init_app, max_iterations, eps, nullptr, std::vector())); + else if (solver_num == 2) + dirichlet_task.set_solver(std::make_unique(init_app, max_iterations, eps, nullptr, std::vector(), Mmin, Mmax)); + else + dirichlet_task.set_solver(std::make_unique(init_app, max_iterations, eps, nullptr, std::vector())); + + auto start = std::chrono::high_resolution_clock::now(); + auto solution = dirichlet_task.solve(); + auto end = std::chrono::high_resolution_clock::now(); + + int64_t duration = std::chrono::duration_cast(end - start).count(); + std::cout << "Overall time first: " << duration / 1000000 << "s\n"; + + std::vector lines; + std::ifstream file("../files/Solver_results.txt"); + if (file.is_open()) + { + std::string line; + while (getline(file, line)) + lines.push_back(line); + file.close(); + } + + n *= 2; + m *= 2; + + std::vector init_app_2((n - 1) * (m - 1), 0.0); + numcpp::DirichletProblemSolver dirichlet_task_2; + + dirichlet_task_2.set_fraction(m, n); + dirichlet_task_2.set_corners(corners); + dirichlet_task_2.set_f(f); + dirichlet_task_2.set_task_num(1); + dirichlet_task_2.set_boundary_conditions({mu1, mu2, mu3, mu4}); + + h = (corners[2] - corners[0]) / n; + k = (corners[3] - corners[1]) / m; + 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); + 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); + + if(solver_num == 0) + dirichlet_task_2.set_solver(std::make_unique + (init_app_2, max_iterations, eps, nullptr, std::vector(), f, mu1, mu2, mu3, mu4, n, m, corners, omega)); + else if(solver_num == 1) + dirichlet_task_2.set_solver(std::make_unique(init_app_2, max_iterations, eps, nullptr, std::vector())); + else if (solver_num == 2) + dirichlet_task_2.set_solver(std::make_unique(init_app_2, max_iterations, eps, nullptr, std::vector(), Mmin, Mmax)); + else + dirichlet_task_2.set_solver(std::make_unique(init_app_2, max_iterations, eps, nullptr, std::vector())); + + auto start_2 = std::chrono::high_resolution_clock::now(); + auto solution_2 = dirichlet_task_2.solve(); + auto end_2 = std::chrono::high_resolution_clock::now(); + + auto duration_2 = std::chrono::duration_cast(end_2 - start_2).count(); + std::cout << "Overall time second: " << duration_2 / 1000000 << "s\n"; + + FP error = 0; + FP x_max = 0; + FP y_max = 0; + + n /= 2; + m /= 2; + + std::ofstream fapp("../files/Approximation.txt"); + std::ofstream fcurr("../files/Correct.txt"); + + for(size_t j = m; j > 0; j--) + { + for(size_t i = 0; i <= n; i++) + { + FP y = corners[1] + i * k*2; + FP x = corners[0] + j * h*2; + FP dot_solution = solution[i][j]; + FP dot_correct = solution_2[i*2][j*2]; + if(std::abs(dot_correct - dot_solution) > error) + { + error = std::abs(dot_correct - dot_solution); + x_max = x; + y_max = y; + } + fapp << dot_solution <<' '; + fcurr << dot_correct << ' '; + } + fapp << '\n'; + fcurr << '\n'; + } + fapp.close(); + fcurr.close(); + + std::ofstream fout_res("../files/Error.txt"); + fout_res << "Задача решена с точностью " << error << "\nМаксимальное отклонение 'точного' и численного решений в точке x = " << x_max << ", y = " << y_max << '\n'; + fout_res.close(); + + std::vector lines_2; + std::ifstream file_2("../files/Solver_results.txt"); + if (file_2.is_open()) + { + std::string line; + while (getline(file_2, line)) + lines_2.push_back(line); + file_2.close(); + } + + std::ofstream fout("../files/Solver_results.txt"); + fout << lines[0] << '\n' << lines[1] << '\n' << lines[2]; + fout.close(); + + std::ofstream fout_2("../files/Solver_results_2N.txt"); + fout_2 << lines_2[0] << '\n' << lines_2[1] << '\n' << lines_2[2]; + fout_2.close(); +} diff --git a/Poisson_Equation/sln/CmakeLists.txt b/Poisson_Equation/sln/CMakeLists.txt similarity index 93% rename from Poisson_Equation/sln/CmakeLists.txt rename to Poisson_Equation/sln/CMakeLists.txt index e29d73a..8765a69 100644 --- a/Poisson_Equation/sln/CmakeLists.txt +++ b/Poisson_Equation/sln/CMakeLists.txt @@ -16,4 +16,4 @@ set(CMAKE_C_FLAGS "${CMAKE_C_FLAGS} -mavx2 -mfma") set(CMAKE_CXX_FLAGS "${CMAKE_CXX_FLAGS} -mavx2 -mfma") # BUILD -add_subdirectory(Dirichlet_Problem) \ No newline at end of file +add_subdirectory(Dirichlet_Problem) diff --git a/Poisson_Equation/sln/Dirichlet_Problem/CmakeLists.txt b/Poisson_Equation/sln/Dirichlet_Problem/CMakeLists.txt similarity index 100% rename from Poisson_Equation/sln/Dirichlet_Problem/CmakeLists.txt rename to Poisson_Equation/sln/Dirichlet_Problem/CMakeLists.txt diff --git a/Poisson_Equation/src/Interface_table.py b/Poisson_Equation/src/Interface_table.py new file mode 100644 index 0000000..2477835 --- /dev/null +++ b/Poisson_Equation/src/Interface_table.py @@ -0,0 +1,405 @@ +import ctypes +import sys +import numpy as np +from PyQt5 import QtWidgets +from PyQt5.QtWidgets import * +from PyQt5.QtCore import Qt + +lib = ctypes.CDLL('./lib_dirichlet.so') +lib.solve_test_task.argtypes = [ctypes.c_uint8, ctypes.c_uint32, ctypes.c_uint32, ctypes.c_uint32, ctypes.c_double, ctypes.c_double] +lib.solve_test_custom_task.argtypes = [ctypes.c_uint8, ctypes.c_uint32, ctypes.c_uint32, ctypes.c_uint32, ctypes.c_double, ctypes.c_double] +lib.solve_main_task.argtypes = [ctypes.c_uint8, ctypes.c_uint32, ctypes.c_uint32, ctypes.c_uint32, ctypes.c_double, ctypes.c_double] + +def test_task(solver_num:int, n: int, m:int, max_iterations:int, eps:float, omega:float): + lib.solve_test_task(solver_num, n, m, max_iterations, eps, omega) + +def test_custom_task(solver_num:int, n: int, m:int, max_iterations:int, eps:float, omega:float): + lib.solve_test_custom_task(solver_num, n, m, max_iterations, eps, omega) + +def main_task(solver_num:int, n: int, m:int, max_iterations:int, eps:float, omega:float): + lib.solve_main_task(solver_num, n, m, max_iterations, eps, omega) + + + +class Window(QMainWindow): + + def __init__(self): + super().__init__() + + self.setWindowTitle("Справка") + self.setGeometry(100, 100, 1040, 1000) + + self.report = QtWidgets.QLabel(self) + self.report.setGeometry(360, 10, 670, 260) + self.report.setStyleSheet('background-color: white;') + + + #making tables + self.table1 = QTableWidget(self) + self.table1.setGeometry(360, 300, 670, 335) + self.table1.setColumnCount(0) + self.table1.setRowCount(0) + + self.table2 = QTableWidget(self) + self.table2.setGeometry(360, 660, 670, 335) + self.table2.setColumnCount(0) + self.table2.setRowCount(0) + + label1 = QtWidgets.QLabel(self) + label1.setText("Аппроксимация") + label1.setGeometry(640, 270, 100, 30) + + label2 = QtWidgets.QLabel(self) + label2.setText("Точное решение") + label2.setGeometry(640, 630, 120, 30) + + + + #list of methods + text1 = QtWidgets.QLabel(self) + text1.setText("Метод") + text1.setGeometry(85, 0, 130, 40) + self.combo_box = QComboBox(self) + self.combo_box.addItem("Метод верхней релаксации") + self.combo_box.addItem("Метод минимальных невязок") + self.combo_box.addItem("Метод Чебышёва") + self.combo_box.addItem("Метод сопряженных градиентов") + self.combo_box.setGeometry(10, 40, 210, 40) + + + text2 = QtWidgets.QLabel(self) + text2.setText("Задача") + text2.setGeometry(85, 120, 100, 40) + self.combo_box1 = QComboBox(self) + self.combo_box1.addItem("Тестовая для стандартной области") + self.combo_box1.addItem("Основная для стандартной области") + self.combo_box1.addItem("Тестовая для нестандартной области") + self.combo_box1.setGeometry(10, 165, 270, 35) + + + text_omega = QtWidgets.QLabel(self) + text_omega.setText("Омега (w) = ") + text_omega.setGeometry(10, 90, 80, 30) + self.omega = QLineEdit(self) + self.omega.setGeometry(90, 90, 50, 30) + self.omega.setText('1.0') + + set1 = QtWidgets.QLabel(self) + set1.setText("Δu(x, y) = ") + set1.setGeometry(10, 250, 80, 30) + self.f = QtWidgets.QLabel(self) + self.f.setGeometry(90, 250, 230, 30) + self.f.setStyleSheet('background-color: white;') + + set2 = QtWidgets.QLabel(self) + set2.setText(" при x ∈ ( -1, 1), y ∈ ( -1, 1);") + set2.setGeometry(10, 280, 200, 30) + + set3 = QtWidgets.QLabel(self) + set3.setText(" u ( -1, y) = ") + set3.setGeometry(10, 310, 80, 30) + self.mu1 = QtWidgets.QLabel(self) + self.mu1.setGeometry(90, 310, 80, 30) + self.mu1.setStyleSheet('background-color: white;') + + + set4 = QtWidgets.QLabel(self) + set4.setText(" u ( 1, y) = ") + set4.setGeometry(170, 310, 80, 30) + self.mu2 = QtWidgets.QLabel(self) + self.mu2.setGeometry(240, 310, 80, 30) + self.mu2.setStyleSheet('background-color: white;') + + + self.set5 = QtWidgets.QLabel(self) + self.set5.setText(" при y ∈ [ -1, 1]; при y ∈ [ -1, 1];") + self.set5.setGeometry(40, 340, 350, 30) + + + set6 = QtWidgets.QLabel(self) + set6.setText(" u ( x, -1) = ") + set6.setGeometry(10, 370, 80, 30) + self.mu3 = QtWidgets.QLabel(self) + self.mu3.setGeometry(90, 370, 80, 30) + self.mu3.setStyleSheet('background-color: white;') + + + set7 = QtWidgets.QLabel(self) + set7.setText(" u ( x, 1) = ") + set7.setGeometry(170, 370, 80, 30) + self.mu4 = QtWidgets.QLabel(self) + self.mu4.setGeometry(240, 370, 80, 30) + self.mu4.setStyleSheet('background-color: white;') + + + self.set8 = QtWidgets.QLabel(self) + self.set8.setText(" при x ∈ [ -1, 1]; при x ∈ [ -1, 1];") + self.set8.setGeometry(40, 400, 350, 30) + + set_line = QtWidgets.QLabel(self) + set_line.setText("__________________________________________________________________________") + set_line.setGeometry(30, 430, 300, 20) + + + set_line = QtWidgets.QLabel(self) + set_line.setText("Доп. ГУ для нестандартной сетки") + set_line.setGeometry(50, 450, 250, 20) + + + set9 = QtWidgets.QLabel(self) + set9.setText(" u ( x, 0) = ") + set9.setGeometry(10, 470, 80, 30) + self.mu5 = QtWidgets.QLabel(self) + self.mu5.setGeometry(90, 470, 80, 30) + self.mu5.setStyleSheet('background-color: white;') + + + set10 = QtWidgets.QLabel(self) + set10.setText(" u ( 0, y) = ") + set10.setGeometry(170, 470, 80, 30) + self.mu6 = QtWidgets.QLabel(self) + self.mu6.setGeometry(240, 470, 80, 30) + self.mu6.setStyleSheet('background-color: white;') + + + set11 = QtWidgets.QLabel(self) + set11.setText(" при x ∈ [ -1, 0];") + set11.setGeometry(40, 500, 120, 30) + + set11 = QtWidgets.QLabel(self) + set11.setText(" при y ∈ [ -1, 0];") + set11.setGeometry(170, 500, 120, 30) + + set_line1 = QtWidgets.QLabel(self) + set_line1.setText("__________________________________________________________________________") + set_line1.setGeometry(30, 530, 300, 20) + + setn = QtWidgets.QLabel(self) + setn.setText("n = ") + setn.setGeometry(10, 560, 30, 30) + setm = QtWidgets.QLabel(self) + setm.setText("m = ") + setm.setGeometry(150, 560, 30, 30) + self.line_editn = QLineEdit(self) + self.line_editn.setGeometry(40, 560, 100, 30) + self.line_editn.setText('8') + self.line_editm = QLineEdit(self) + self.line_editm.setGeometry(180, 560, 100, 30) + self.line_editm.setText('8') + + + setit = QtWidgets.QLabel(self) + setit.setText("Максимальное количество итераций ") + setit.setGeometry(10, 600, 240, 30) + self.line_editit = QLineEdit(self) + self.line_editit.setGeometry(250, 600, 100, 30) + self.line_editit.setText('1000000') + + seteps = QtWidgets.QLabel(self) + seteps.setText("Требуемая точность решения СЛАУ") + seteps.setGeometry(10, 635, 230, 30) + self.line_editeps = QLineEdit(self) + self.line_editeps.setGeometry(240, 635, 110, 30) + self.line_editeps.setText('0.00000001') + + + #button set task + self.button1 = QtWidgets.QPushButton(self) + self.button1.setGeometry(15, 210, 200, 30) + self.button1.setText("Задать условия задачи") + self.button1.clicked.connect(self.Set_task) + + + #button solve task + self.button2 = QtWidgets.QPushButton(self) + self.button2.setGeometry(30, 680, 200, 30) + self.button2.setText("Аппроксимировать") + self.button2.clicked.connect(self.Solve_task) + + def Set_task(self): + if(self.combo_box1.currentIndex() == 0): + self.f.setText("4exp(1 - x^2 - y^2)(x^2 + y^2 - 1)") + self.mu1.setText("exp(- y^2)") + self.mu2.setText("exp(- y^2)") + self.mu3.setText("exp(- x^2)") + self.mu4.setText("exp(- x^2)") + self.mu5.setText("") + self.mu6.setText("") + self.set5.setText(" при y ∈ [ -1, 1]; при y ∈ [ -1, 1];") + self.set8.setText(" при x ∈ [ -1, 1]; при x ∈ [ -1, 1];") + + + elif(self.combo_box1.currentIndex() == 1): + self.f.setText("- |(sin(πxy))^3|") + self.mu1.setText("1- y^2") + self.mu2.setText("1 - y^2") + self.mu3.setText("1 - x^2") + self.mu4.setText("1 - x^2") + self.mu5.setText("") + self.mu6.setText("") + self.set5.setText(" при y ∈ [ -1, 1]; при y ∈ [ -1, 1];") + self.set8.setText(" при x ∈ [ -1, 1]; при x ∈ [ -1, 1];") + else: + self.f.setText("4exp(1 - x^2 - y^2)(x^2 + y^2 - 1)") + self.mu1.setText("exp(- y^2)") + self.mu2.setText("exp(- y^2)") + self.mu3.setText("exp(- x^2)") + self.mu4.setText("exp(- x^2)") + self.mu5.setText("exp(1 - x^2)") + self.mu6.setText("exp(1 - y^2)") + self.set5.setText(" при y ∈ [ 0, 1]; при y ∈ [ -1, 1];") + self.set8.setText(" при x ∈ [ 0, 1]; при x ∈ [ -1, 1];") + + + def Solve_task(self): + n = (int)(self.line_editn.text()) + m = (int)(self.line_editm.text()) + omega = (float)(self.omega.text()) + max_iterations = (int)(self.line_editit.text()) + eps = (float)(self.line_editeps.text()) + task_num = self.combo_box1.currentIndex() + solver_num = (int)(self.combo_box.currentIndex()) + solver_str = " " + if(task_num == 0): + solver_str += "Тестовая задача на стандартной сетке решена " + elif(task_num == 1): + solver_str += "Основная задача на стандартной сетке решена " + else: + solver_str += "Тестовая задача на нестандартной сетке решена " + + if(solver_num == 0): + solver_str += "методом верхней релаксации \nс параметром w = "+ self.omega.text() + elif(solver_num == 1): + solver_str += "методом минимальных невязок" + elif(solver_num == 2): + solver_str += "методом Чебышева" + else: + solver_str += "методом сопряженных градиентов" + solver_str += "\n с нулевым начальным приближением " + + if(task_num == 0): + test_task(solver_num, n, m, max_iterations, eps, omega) + + error = open("../files/Error.txt", 'r') + solver_results = open("../files/Solver_results.txt", 'r') + results = [row.strip() for row in solver_results] + + self.report.setText("Для решения тестовой задачи использованы сетка с числом разбиений\n по х : n = " + (str)(n)+ + ", и числом разбиений по y: m = "+(str)(m)+",\n " + solver_str + ", применены критерии остановки \n по точности решения СЛАУ eps(мет) = "+ + (str)(eps)+" и по числу итераций N(max) = "+ (str)(max_iterations)+".\n \n " + "На решение схемы (СЛАУ) затрачено " + (str)(results[2]) + + " итераций и достигнута точность " + (str)(results[0]) + "\nСхема(СЛАУ) решена с невязкой по норме Чебышёва " + (str)(results[1]) + "\n" + error.read()) + error.close() + solver_results.close() + + + elif(task_num == 1): + main_task(solver_num, n, m, max_iterations, eps, omega) + + curr_string = "\nДля контроля точности использована сетка(2N) \n с числом разбиений по x : n = " + curr_string += (str)(n*2) + curr_string += ", и числом разбиений по y: m = " + curr_string += (str)(m*2) + curr_string += ". \nКритерии остановки метода остаются такими же." + + error = open("../files/Error.txt", 'r') + solver_results = open("../files/Solver_results.txt", 'r') + results = [row.strip() for row in solver_results] + + solver_results_2 = open("../files/Solver_results_2N.txt", 'r') + results_2 = [row.strip() for row in solver_results_2] + + curr_string += "На решение СЛАУ(2N) затрачено " + curr_string += (str)(results_2[2]) + curr_string += " итераций\n и достигнута точность " + curr_string += (str)(results_2[0]) + curr_string += ".СЛАУ(2N) решена с невязкой по норме Чебышёва " + curr_string += (str)(results_2[1]) + + self.report.setText("Для решения основной задачи использованы сетка с числом разбиений\n по х : n = " + (str)(n)+ + ", и числом разбиений по y: m = "+(str)(m)+",\n " + solver_str + ", применены критерии остановки \n по точности решения СЛАУ eps(мет) = "+ + (str)(eps)+" и по числу итераций N(max) = "+ (str)(max_iterations)+".\n \n " + "На решение схемы (СЛАУ) затрачено " + (str)(results[2]) + + " итераций и достигнута точность " + (str)(results[0]) + "\nСхема(СЛАУ) решена с невязкой по норме Чебышёва " + (str)(results[1]) + "\n" + curr_string + "\n\n" + error.read()) + error.close() + solver_results.close() + + else: + test_custom_task(solver_num, n, m, max_iterations, eps, omega) + + error = open("../files/Error.txt", 'r') + solver_results = open("../files/Solver_results.txt", 'r') + results = [row.strip() for row in solver_results] + + self.report.setText("Для решения тестовой задачи использованы сетка с числом разбиений\n по х : n = " + (str)(n)+ + ", и числом разбиений по y: m = "+(str)(m)+",\n " + solver_str + ", применены критерии остановки \n по точности решения СЛАУ eps(мет) = "+ + (str)(eps)+" и по числу итераций N(max) = "+ (str)(max_iterations)+".\n \n " + "На решение схемы (СЛАУ) затрачено " + (str)(results[2]) + + " итераций и достигнута точность " + (str)(results[0]) + "\nСхема(СЛАУ) решена с невязкой по норме Чебышёва " + (str)(results[1]) + "\n" + error.read()) + error.close() + solver_results.close() + self.load_data_to_table1("../files/Approximation.txt") + self.load_data_to_table2("../files/Correct.txt") + + def load_data_to_table1(self, filename): + self.table1.setRowCount(0) + self.table1.setColumnCount(0) + try: + with open(filename, 'r') as file: + lines = file.readlines() + + # Определяем количество строк и столбцов + num_rows = len(lines) + num_cols = max(len(line.split()) for line in lines) + + # Устанавливаем количество строк и столбцов в таблице + self.table1.setRowCount(num_rows) + self.table1.setColumnCount(num_cols) + + # Заполняем таблицу данными + for row_index, line in enumerate(lines): + values = line.split() + for col_index, value in enumerate(values): + item = QTableWidgetItem(value) + self.table1.setItem(row_index, col_index, item) + except FileNotFoundError: + print("Файл не найден") + except Exception as e: + print(f"Ошибка при чтении файла: {e}") + + def load_data_to_table2(self, filename): + self.table2.setRowCount(0) + self.table2.setColumnCount(0) + try: + with open(filename, 'r') as file: + lines = file.readlines() + + # Определяем количество строк и столбцов + num_rows = len(lines) + num_cols = max(len(line.split()) for line in lines) + + # Устанавливаем количество строк и столбцов в таблице + self.table2.setRowCount(num_rows) + self.table2.setColumnCount(num_cols) + + # Заполняем таблицу данными + for row_index, line in enumerate(lines): + values = line.split() + for col_index, value in enumerate(values): + item = QTableWidgetItem(value) + self.table2.setItem(row_index, col_index, item) + except FileNotFoundError: + print("Файл не найден") + except Exception as e: + print(f"Ошибка при чтении файла: {e}") + + + + + + + + + +if __name__ == '__main__': + app = QApplication(sys.argv) + window = Window() + window.show() + sys.exit(app.exec_()) \ No newline at end of file diff --git a/Poisson_Equation/src/compile_lib.sh b/Poisson_Equation/src/compile_lib.sh new file mode 100644 index 0000000..0dbbafc --- /dev/null +++ b/Poisson_Equation/src/compile_lib.sh @@ -0,0 +1 @@ +g++ -fPIC -shared -mavx2 -mfma -fopenmp -std=c++17 -O3 utility.cpp -o lib_dirichlet.so diff --git a/Poisson_Equation/src/lib_dirichlet.so b/Poisson_Equation/src/lib_dirichlet.so new file mode 100644 index 0000000..3b02f49 Binary files /dev/null and b/Poisson_Equation/src/lib_dirichlet.so differ diff --git a/Poisson_Equation/src/main.cpp b/Poisson_Equation/src/main.cpp index 07aa420..fc662d0 100644 --- a/Poisson_Equation/src/main.cpp +++ b/Poisson_Equation/src/main.cpp @@ -1,5 +1,17 @@ -#include "Dirichlet_Problem.hpp" -#include "Solvers.hpp" +#include +#include +#include +#include +#include +#include +#include +#include + +#include "../headers/Dirichlet_Problem.hpp" +#include "../headers/Solvers.hpp" + +#define FP double + int main() { return 0; diff --git a/Poisson_Equation/src/test b/Poisson_Equation/src/test new file mode 100755 index 0000000..3c5cd11 Binary files /dev/null and b/Poisson_Equation/src/test differ diff --git a/Poisson_Equation/src/test.cpp b/Poisson_Equation/src/test.cpp index 18b37ce..d092efc 100644 --- a/Poisson_Equation/src/test.cpp +++ b/Poisson_Equation/src/test.cpp @@ -6,9 +6,13 @@ #include #include #include +#include +#include +#include #include "../headers/Dirichlet_Problem.hpp" #include "../headers/Solvers.hpp" +#include "../headers/tasks.hpp" #define FP double @@ -25,105 +29,116 @@ std::pair>> estimate_time(numcpp::DirichletP return {duration, solution}; } -void test_TopRelaxation(){ - size_t n = 2048; - size_t m = 2048; - std::cout << "n = " << n << " m = " << m << '\n'; - - FP omega = 1.95; - std::array corners = {-1.0, -1.0, 1.0, 1.0}; - std::vector init_app((n - 1) * (m - 1), 0.0); - - 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)); }; +void report(size_t task_num, size_t n, size_t m, size_t max_iterations, FP eps) +{ + std::vector results; + std::vector results_err; + std::string line = ""; + std::string curr_string = ""; + std::string start_string = ""; - 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}); + std::ifstream solver_results("../files/Solver_results.txt"); + while (std::getline(solver_results, line)) + results.push_back(line); + solver_results.close(); - //dirichlet_task.set_solver(std::make_unique(init_app, 10000, 0.00001, nullptr, std::vector(), omega)); - dirichlet_task.set_solver(std::make_unique - (init_app, 1000000000, 0.0000000000001, nullptr, std::vector(), f, mu1, mu2, mu3, mu4, n, m, corners, omega)); + std::ifstream error("../files/Error.txt"); + while (std::getline(error, line)) + results_err.push_back(line); + error.close(); - auto [duration, solution] = estimate_time(dirichlet_task); -} -void test_task(size_t m, size_t n, std::unique_ptr LS_solver) -{ - std::cout << "n = " << n << " m = " << m << '\n'; - size_t size = (n - 1) * (m - 1); - std::array corners = {-1.0, -1.0, 1.0, 1.0}; - - std::vector initial_approximation(size, 0.0); - LS_solver->set_initial_approximation(initial_approximation); - - 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)); }; + if(task_num == 0) + { + start_string = "Для решения тестовой задачи использованы сетка "; + } + else if(task_num == 1) + { + start_string = "Для решения основной задачи использованы сетка "; + + curr_string = "\nДля контроля точности использована сетка(2N) \n с числом разбиений по x : n = "; + curr_string += std::to_string(n * 2); + curr_string += ", и числом разбиений по y: m = "; + curr_string += std::to_string(m * 2); + curr_string += ". \nКритерии остановки метода остаются такими же."; + + std::ifstream solver_results_2("../files/Solver_results_2N.txt"); + std::vector results_2; + while (std::getline(solver_results_2, line)) + results_2.push_back(line); + solver_results_2.close(); + + + curr_string += "На решение СЛАУ(2N) затрачено "; + curr_string += results_2[2]; + curr_string += " итераций\n и достигнута точность "; + curr_string += results_2[0]; + curr_string += ".СЛАУ(2N) решена с невязкой по норме Чебышёва "; + curr_string += results_2[1]; + } + else + { + start_string = "Для решения тестовой задачи использованы сетка (нестандартная сетка) "; + } + + std::cout << "\n\n\n\n" << start_string << " с числом разбиений\n по х : n = " << n + << ", и числом разбиений по y: m = " << m << ",\n " + << "применены критерии остановки \n по точности решения СЛАУ eps(мет) = " << eps + << " и по числу итераций N(max) = " << max_iterations << ".\n \n " + << "На решение схемы (СЛАУ) затрачено " << results[2] + << " итераций и достигнута точность " << results[0] + << "\nСхема(СЛАУ) решена с невязкой по норме Чебышёва " << results[1] << "\n" << results_err[0] << "\n" << results_err[1] << "\n\n" << curr_string << "\n\n\n"; - numcpp::DirichletProblemSolver solver; - solver.set_fraction(m, n); - solver.set_corners(corners); - solver.set_u(u); - solver.set_f(f); - solver.set_boundary_conditions({mu1, mu2, mu3, mu4}); - solver.set_solver(std::move(LS_solver)); - - auto [duration, solution] = estimate_time(solver); } -void test_task_custom_grid(size_t m, size_t n, std::unique_ptr LS_solver) -{ - size_t size = (n / 2 - 1) * (m - 1) + (n / 2) * (m / 2 - 1); - std::cout << "n = " << n << " m = " << m << '\n'; - - std::vector initial_approximation(size, 0.0); - LS_solver->set_initial_approximation(initial_approximation); - - std::array corners = {-1.0, -1.0, 1.0, 1.0}; - - 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(1.0 - pow(y, 2)); }; - auto mu3 = [](double y) { return exp(-pow(y, 2)); }; - auto mu4 = [](double x) { return exp(-pow(x, 2)); }; - auto mu5 = [](double x) { return exp(1.0 - pow(x, 2)); }; - auto mu6 = [](double x) { return exp(-pow(x, 2)); }; - - numcpp::DirichletProblemSolver solver; - std::array, 6> arr{mu1, mu2, mu3, mu4, mu5, mu6}; - solver.set_solver(std::move(LS_solver)); - solver.set_boundary_conditions_for_r_shaped_grid(arr); - solver.set_fraction(n, m); - solver.set_corners(corners); - solver.set_u(u); - solver.set_f(f); - - auto [duration, solution] = estimate_time(solver); -} +// for test we use functions from tasks.hpp int main() { - size_t m = 1024; - size_t n = 1024; - - auto LS_solver = std::make_unique(std::vector(), 1000000, 0.00001, nullptr, std::vector()); - - test_task_custom_grid(m, n, std::move(LS_solver)); -} + size_t m = 500; + size_t n = 500; + // solver_num chooses a method for solving a system of linear equations + // 0 = TopRelaxation + // 1 = MinRes + // 2 = ChebyshevIteration + // 3 = ConGrad + //size_t solver_num = 3; + size_t max_iterations = 1000000; + FP eps = 0.00000001; + FP omega = 1.99; + // task_num + // 0 = test + // 1 = main + // 2 = custom test + size_t task_num = 0; + + std::cout << "========== TopRelaxation ==========" << std::endl; + test_task(0, n, m, max_iterations, eps, omega); // Regular grid + report(task_num, n, m, max_iterations, eps); + std::cout << std::endl; + + std::cout << "========== MinRes ==========" << std::endl; + test_task(1, n, m, max_iterations, eps, omega); // Regular grid + report(task_num, n, m, max_iterations, eps); + std::cout << std::endl; + + std::cout << "========== ChebyshevIteration ==========" << std::endl; + test_task(2, n, m, max_iterations, eps, omega); // Regular grid + report(task_num, n, m, max_iterations, eps); + std::cout << std::endl; + + std::cout << "========== ConGrad ==========" << std::endl; + test_task(3, n, m, max_iterations, eps, omega); // Regular grid + report(task_num, n, m, max_iterations, eps); + std::cout << std::endl; + + //main_task(solver_num, n, m, max_iterations, eps, omega); // Regular grid + // test_custom_task(solver_num, n, m, max_iterations, eps, omega); //ReversedR grid + + // the values in the grid nodes can be found in the files + // for approximation in file "../files/Approximation.txt" + // for correct solution in file "../files/Correct.txt" +} \ No newline at end of file diff --git a/Poisson_Equation/src/utility.cpp b/Poisson_Equation/src/utility.cpp new file mode 100644 index 0000000..06eee14 --- /dev/null +++ b/Poisson_Equation/src/utility.cpp @@ -0,0 +1,20 @@ +#include "../headers/tasks.hpp" + +#define FP double +static uint8_t solver_num; +static uint32_t n, m, max_iterations; +static FP eps, omega; + +extern "C" void solve_test_task(uint8_t solver_num, uint32_t n, uint32_t m, uint32_t max_iterations, FP eps, FP omega) { + test_task(solver_num, n, m, max_iterations, eps, omega); +} + + +extern "C" void solve_test_custom_task(uint8_t solver_num, uint32_t n, uint32_t m, uint32_t max_iterations, FP eps, FP omega) { + test_custom_task(solver_num, n, m, max_iterations, eps, omega); +} + + +extern "C" void solve_main_task(uint8_t solver_num, uint32_t n, uint32_t m, uint32_t max_iterations, FP eps, FP omega) { + main_task(solver_num, n, m, max_iterations, eps, omega); +} \ No newline at end of file