From 7ec31ae16963aeeb544aaba050b1e7df82331a54 Mon Sep 17 00:00:00 2001 From: Allan K <73891596+CodingAllan@users.noreply.github.com> Date: Mon, 17 Jun 2024 10:02:27 +0200 Subject: [PATCH] 33 prolongation and restriction tests (#34) MIME-Version: 1.0 Content-Type: text/plain; charset=UTF-8 Content-Transfer-Encoding: 8bit Add tests for prolongation and restriction operators. Co-authored-by: Martin J. Kühn <62713180+mknaranja@users.noreply.github.com> Co-authored-by: Martin Joachim Kühn --- tests/CMakeLists.txt | 5 +- tests/mockgrid.cpp | 67 ++++++++ tests/mockgrid.h | 11 ++ tests/test_prolongation.cpp | 299 ++++++++++++++++++++++++++++++++++++ tests/test_restriction.cpp | 295 +++++++++++++++++++++++++++++++++++ 5 files changed, 675 insertions(+), 2 deletions(-) create mode 100644 tests/mockgrid.cpp create mode 100644 tests/mockgrid.h create mode 100644 tests/test_prolongation.cpp create mode 100644 tests/test_restriction.cpp diff --git a/tests/CMakeLists.txt b/tests/CMakeLists.txt index 7fe03ba5..90fc6b85 100644 --- a/tests/CMakeLists.txt +++ b/tests/CMakeLists.txt @@ -1,6 +1,6 @@ enable_testing() -set(TESTSOURCES comparetest.cpp) +set(TESTSOURCES comparetest.cpp test_prolongation.cpp test_restriction.cpp mockgrid.cpp) add_executable( gmgpolar_tests @@ -12,6 +12,7 @@ target_link_libraries( GMGPolar ) -target_include_directories(gmgpolar_tests PRIVATE ../include ../include/test_cases) + +target_include_directories(gmgpolar_tests PRIVATE . ../include ../include/test_cases) include(GoogleTest) gtest_discover_tests(gmgpolar_tests) \ No newline at end of file diff --git a/tests/mockgrid.cpp b/tests/mockgrid.cpp new file mode 100644 index 00000000..4efe7a81 --- /dev/null +++ b/tests/mockgrid.cpp @@ -0,0 +1,67 @@ +#include "mockgrid.h" +void create_grid(gmgpolar& test_p) +{ + time_t seed = 5000; + std::default_random_engine gen(seed); //deterministic seed to reproduce + std::uniform_real_distribution dis(gyro::dcntl[Param::R0], gyro::dcntl[Param::R]); + std::uniform_real_distribution theta_distribution(0, 2 * PI); + + std::cout << "creating randomized grid with std::default_random_engine seed " + std::to_string(seed) << std::endl; + level* new_level = new level(0); + new_level->nr = pow(2, gyro::icntl[Param::nr_exp]); + new_level->r = std::vector(new_level->nr + 1); + + std::vector rands(new_level->nr - 1); + for (int j = 0; j < new_level->nr - 1; ++j) { + rands[j] = dis(gen); + } + std::sort(rands.begin(), rands.end()); + new_level->r[0] = gyro::dcntl[Param::R0]; + for (int i = 1; i < new_level->nr; i++) { + new_level->r[i] = rands[i - 1]; //random grid on coarser level + } + new_level->r[new_level->nr] = gyro::dcntl[Param::R]; + new_level->nr++; + int ntmp = pow(2, ceil(log2(new_level->nr))); + new_level->ntheta = gyro::icntl[Param::periodic] ? ntmp : ntmp + 1; + + new_level->theta = std::vector(new_level->ntheta); + + std::vector randst(ntmp - 1); + for (int k = 0; k < ntmp - 1; ++k) { + randst[k] = theta_distribution(gen); + } + std::sort(randst.begin(), randst.end()); + new_level->theta[0] = 0; + for (int i = 1; i < ntmp; i++) { + new_level->theta[i] = randst[i - 1]; + } + if (!gyro::icntl[Param::periodic]) { + new_level->theta[ntmp] = 2 * PI; + } + + new_level->ntheta_int = gyro::icntl[Param::periodic] ? new_level->ntheta : new_level->ntheta - 1; + new_level->nr_int = new_level->nr - 1; + + new_level->thetaplus = std::vector(new_level->ntheta_int); + for (int k = 0; k < new_level->ntheta_int - 1; k++) { + new_level->thetaplus[k] = new_level->theta[k + 1] - new_level->theta[k]; + } + new_level->thetaplus[new_level->ntheta_int - 1] = 2 * PI - new_level->theta[new_level->ntheta_int - 1]; + + new_level->hplus = std::vector(new_level->nr_int); + for (int k = 0; k < new_level->nr_int; k++) { + new_level->hplus[k] = new_level->r[k + 1] - new_level->r[k]; + } + + level* coarser_level = new level(1); + coarser_level->nr = pow(2, gyro::icntl[Param::nr_exp] - 1) + 1; + ntmp = pow(2, ceil(log2(coarser_level->nr))); + coarser_level->ntheta = gyro::icntl[Param::periodic] ? ntmp : ntmp + 1; + + coarser_level->ntheta_int = gyro::icntl[Param::periodic] ? coarser_level->ntheta : coarser_level->ntheta - 1; + coarser_level->nr_int = coarser_level->nr - 1; + test_p.v_level.push_back(new_level); + test_p.v_level.push_back(coarser_level); + test_p.levels_orig = 2; +} \ No newline at end of file diff --git a/tests/mockgrid.h b/tests/mockgrid.h new file mode 100644 index 00000000..0d000a14 --- /dev/null +++ b/tests/mockgrid.h @@ -0,0 +1,11 @@ +#include "gmgpolar.h" +#include +#include +#include + +#ifndef MOCKGRID_H + #define MOCKGRID_H + +void create_grid(gmgpolar& test_p); + +#endif \ No newline at end of file diff --git a/tests/test_prolongation.cpp b/tests/test_prolongation.cpp new file mode 100644 index 00000000..a78f4f1a --- /dev/null +++ b/tests/test_prolongation.cpp @@ -0,0 +1,299 @@ +#include +#include +#include "gmgpolar.h" +#include "mockgrid.h" +class test_prolongation + : public ::testing::TestWithParam< + std::tuple> //tuple includes data on grid size and Dirichlet boundary conditions +{ +protected: + void SetUp() override + { + //initialize default parameters. + gyro::init_params(); + gyro::icntl[Param::verbose] = 0; + gyro::dcntl[Param::R0] = 1e-5; + gyro::icntl[Param::periodic] = 1; + gyro::f_grid_r = ""; + gyro::f_grid_theta = ""; + gyro::f_sol_in = ""; + gyro::f_sol_out = ""; + gyro::select_functions_class(gyro::icntl[Param::alpha_coeff], gyro::icntl[Param::beta_coeff], + gyro::icntl[Param::mod_pk], gyro::icntl[Param::prob]); + } +}; + +/*! + * \brief Test the bilinear prolongation operator used in the multigrid cycle coarse-grid correction. + * + * The Test creates an arbitrary grid-function on the coarser level and prolongates it. + * On the fine level we iterate over all nodes and test the result based on whether the node is fine in theta, r or in both. + * + * Parametrized tests are used to test for different grid sizes and with or without Dirichlet boundary conditions. + */ + +TEST_P(test_prolongation, test_bilinear_prolongation) +{ + //we vary the grid size to guarantee that the problem works for all sizes + gyro::icntl[Param::nr_exp] = (int)(std::get<0>(GetParam()) / 3) + 3; + gyro::icntl[Param::ntheta_exp] = (std::get<0>(GetParam()) % 3) + 3; + gyro::icntl[Param::DirBC_Interior] = std::get<1>(GetParam()); + + gmgpolar test_p; + create_grid(test_p); + + level& p_level = *(test_p.v_level[0]); + int ctheta_int = test_p.v_level[1]->ntheta_int; //number of coarse nodes in theta direction + + p_level.m = test_p.v_level[0]->nr * test_p.v_level[0]->ntheta; //fine grid size + p_level.mc = test_p.v_level[1]->nr * test_p.v_level[1]->ntheta; //coarser grid size + + std::vector u_test(p_level.mc); + for (int z = 0; z < p_level.mc; z++) { + u_test[z] = + 1 - z + + pow(PI, + -z * z); //constructing arbitrary grid-function on coarse level to test our prolongation operator with. + } + + std::vector sol = + p_level.apply_prolongation_bi(u_test); //apply prolongation operator on arbitrary grid-function + + for (int j = 0; j < p_level.nr_int + 1; j++) { + for (int i = 0; i < p_level.ntheta_int; i++) { + if (j % 2 == 0 && i % 2 == 0) { //testing coarse node injection + EXPECT_EQ(u_test[(j / 2) * ctheta_int + (i / 2)], sol[j * p_level.ntheta_int + i]) + << "coarse node injection is failing"; + } + else if (i % 2 != 0 && j % 2 == 0) { + // as numbering in angle (theta_i) is odd, we have a fine node with + // coarse neighbors at (r_j, theta_i - k_{i-1}) and (r_j, theta_i + k_i) + + double k_qm1 = p_level.theta[i] - p_level.theta[i - 1]; //calculate k_{q-1} + double k_q = (i < p_level.ntheta_int - 1) ? p_level.theta[i + 1] - p_level.theta[i] + : 2 * PI - p_level.theta[i]; //k_q + + int i1 = (j / 2) * ctheta_int + (i - 1) / 2; //bottom coarse node in theta + int i2 = (i < p_level.ntheta_int - 1) ? i1 + 1 : (j / 2) * ctheta_int; //top coarse node in theta + + double val = (k_q * u_test[i1] + k_qm1 * u_test[i2]); + + EXPECT_NEAR(val, (k_q + k_qm1) * sol[j * p_level.ntheta_int + i], + 1e-12) //compare values as in the paper + << "Bilinear Prolongation fails for Index (r,theta) : (" + std::to_string(j) + "," + + std::to_string(i) + ")"; + ; + } + else if (i % 2 == 0 && j % 2 != 0) { + // as numbering in radius (r_j) is odd, we have a fine node with + // coarse neighbors at (r_j - h_{j-1}, theta_i) and (r_j+h_j, theta_i ) + double h_pm1 = p_level.r[j] - p_level.r[j - 1]; //h_{p-1} + double h_p = p_level.r[j + 1] - p_level.r[j]; //h_p + + double v1 = u_test[(j - 1) / 2 * ctheta_int + (i / 2)]; // left coarse node in r + double v2 = u_test[(j + 1) / 2 * ctheta_int + (i / 2)]; // right coarse node in r + + double val = (h_p * v1 + h_pm1 * v2); + + EXPECT_NEAR(val, (h_p + h_pm1) * sol[j * p_level.ntheta_int + i], 1e-12) //compare values + << "Bilinear Prolongation fails for Index (r,theta) : (" + std::to_string(j) + "," + + std::to_string(i) + ")"; + ; + } + else // both are fine nodes + { + double h_pm1 = p_level.r[j] - p_level.r[j - 1]; + double h_p = p_level.r[j + 1] - p_level.r[j]; + + double k_qm1 = p_level.theta[i] - p_level.theta[i - 1]; + double k_q = + (i < p_level.ntheta_int - 1) ? p_level.theta[i + 1] - p_level.theta[i] : 2 * PI - p_level.theta[i]; + + /*the stencil-corners corresponding to the values (r_j,theta_i)*/ + /*bottom corresponds to lower indices in theta direction. left to lower indices in radius direction*/ + + double bottom_left; + double top_left; + double bottom_right; + double top_right; + ASSERT_NE(p_level.nr_int - 1, 1); + + if (i < p_level.ntheta - 1) { + bottom_left = u_test[((j - 1) / 2) * ctheta_int + (i - 1) / 2]; + top_left = u_test[((j - 1) / 2) * ctheta_int + (i + 1) / 2]; + bottom_right = u_test[((j + 1) / 2) * ctheta_int + (i - 1) / 2]; + top_right = u_test[((j + 1) / 2) * ctheta_int + (i + 1) / 2]; + } + else { + bottom_left = u_test[((j - 1) / 2) * ctheta_int + (i - 1) / 2]; + top_left = u_test[((j - 1) / 2) * ctheta_int]; + bottom_right = u_test[((j + 1) / 2) * ctheta_int + (i - 1) / 2]; + top_right = u_test[((j + 1) / 2) * ctheta_int]; + } + + double val = ((h_p * k_q * bottom_left) + (h_p * k_qm1 * top_left) + (h_pm1 * k_q * bottom_right) + + (h_pm1 * k_qm1 * top_right)); //calculate value as in the paper + + EXPECT_NEAR(val, (h_p + h_pm1) * (k_q + k_qm1) * sol[j * p_level.ntheta_int + i], 1e-12) + << "Bilinear Prolongation fails for Index (r,theta) : (" + std::to_string(j) + "," + + std::to_string(i) + ")"; + } + } + } +} + +/*! + * \brief Test the injection prolongation operator used in the implicit extrapolation step of the multigrid-cycle. + * + * The Test creates an arbitrary grid-function on the coarser level and injects it. + * On the fine level we iterate over all nodes and test the result based on whether the node is fine or not. + * + * Parametrized tests are used to test for different grid sizes and with or without Dirichlet boundary conditions. + */ + +TEST_P(test_prolongation, test_injection_prolongation) +{ + gyro::icntl[Param::nr_exp] = (int)(std::get<0>(GetParam()) / 3) + 3; + gyro::icntl[Param::ntheta_exp] = (std::get<0>(GetParam()) % 3) + 3; + gyro::icntl[Param::DirBC_Interior] = std::get<1>(GetParam()); + + gmgpolar test_p; + + create_grid(test_p); + + level& p_level = *(test_p.v_level[0]); + int ctheta_int = test_p.v_level[1]->ntheta_int; // number of coarse nodes in theta direction + + p_level.m = test_p.v_level[0]->nr * test_p.v_level[0]->ntheta; // fine grid size + p_level.mc = test_p.v_level[1]->nr * test_p.v_level[1]->ntheta; // coarser grid size + + std::vector u_test(p_level.mc); + for (int z = 0; z < p_level.mc; z++) { + u_test[z] = + 1 - z + + pow(PI, + -z * z); //constructing arbitrary grid-function on coarse level to test our prolongation operator with. + } + + std::vector sol = p_level.apply_prolongation_inj(u_test); + + EXPECT_EQ((int)sol.size(), p_level.m); + //testing the embedding. since we embed the coarse gird function, the values of the fine node should be zero + for (int j = 0; j < p_level.nr_int + 1; j++) { + for (int i = 0; i < p_level.ntheta_int; i++) { + if (j % 2 == 0 && i % 2 == 0) { + EXPECT_EQ(u_test[(j / 2) * ctheta_int + (i / 2)], sol[j * p_level.ntheta_int + i]) + << "The Injection value fails at Index (r,theta): (" + std::to_string(j) + "," + std::to_string(i) + + ")"; + } + else { + EXPECT_EQ(sol[j * p_level.ntheta_int + i], 0); //fine nodes set to zero. + } + } + } +} +/*! + * \brief Test the extrapolation prolongation operator used in the implicit extrapolation step of the multigrid-cycle. + * + * The Test creates an arbitrary grid-function on the coarser level and prolongates it. + * On the fine level we iterate over all nodes and test the result based on whether the node is fine in theta, r or in both. + * The value is determined by a Triangulation of the grid on the coarser level. + * + * Parametrized tests are used to test for different grid sizes and with or without Dirichlet boundary conditions. + */ + +TEST_P(test_prolongation, test_extrapolation_prolongation) +{ + gyro::icntl[Param::nr_exp] = (int)(std::get<0>(GetParam()) / 3) + 3; + gyro::icntl[Param::ntheta_exp] = (std::get<0>(GetParam()) % 3) + 3; + gyro::icntl[Param::DirBC_Interior] = std::get<1>(GetParam()); + + gmgpolar test_p; + + create_grid(test_p); + + level& p_level = *(test_p.v_level[0]); + int ctheta_int = test_p.v_level[1]->ntheta_int; // number of coarse nodes in theta direction + + p_level.m = test_p.v_level[0]->nr * test_p.v_level[0]->ntheta; // number of nodes on fine level + p_level.mc = test_p.v_level[1]->nr * test_p.v_level[1]->ntheta; // number of nodes on coarse level + + std::vector u_test(p_level.mc); + for (int z = 0; z < p_level.mc; z++) { + u_test[z] = + 1 - z + + pow(PI, + -z * z); //constructing arbitrary grid-function on coarse level to test our prolongation operator with. + } + std::vector sol = p_level.apply_prolongation_ex(u_test); + + for (int j = 0; j < p_level.nr_int + 1; j++) { + for (int i = 0; i < p_level.ntheta_int; i++) { + if (j % 2 == 0 && i % 2 == 0) { //coarse node injection + EXPECT_EQ(u_test[(j / 2) * ctheta_int + (i / 2)], sol[j * p_level.ntheta_int + i]) + << "coarse node injection is failing"; + } + else if (i % 2 != 0 && j % 2 == 0) { + // as numbering in angle (theta_i) is odd, we have a fine node with + // coarse neighbors at (r_j, theta_i - k_{i-1}) and (r_j, theta_i + k_i) + + int i1 = (j / 2) * ctheta_int + (i - 1) / 2; // bottom coarse node in theta + int i2 = (i < p_level.ntheta_int - 1) ? i1 + 1 : (j / 2) * ctheta_int; // top coarse node in theta + + double val = 0.5 * (u_test[i1] + u_test[i2]); + + EXPECT_NEAR(val, sol[j * p_level.ntheta_int + i], 1e-12) + << "Extrapolated Prolongation fails for Index (r,theta): (" + std::to_string(j) + "," + + std::to_string(i) + ")"; + } + else if (i % 2 == 0 && j % 2 != 0) { + // as numbering in radius (r_j) is odd, we have a fine node with + // coarse neighbors at (r_j - h_{j-1}, theta_i) and (r_j+h_j, theta_i ) + + double v1 = u_test[(j - 1) / 2 * ctheta_int + (i / 2)]; // left coarse node in r + double v2 = u_test[(j + 1) / 2 * ctheta_int + (i / 2)]; // right coarse node in r + + double val = 0.5 * (v1 + v2); + + EXPECT_NEAR(val, sol[j * p_level.ntheta_int + i], 1e-12) + << "Extrapolated Prolongation fails for Index (r,theta): (" + std::to_string(j) + "," + + std::to_string(i) + ")"; + } + else // both are fine nodes + { + + /*in the triangulation we now consider that the fine node is on the hypothenuse of the triangle. + The corresponding coarse nodes are located on the triangles nodes that span the hypothenuse */ + + double top_left; + double bottom_right; + + if (i < p_level.ntheta_int - 1) { + top_left = u_test[((j - 1) / 2) * ctheta_int + (i + 1) / 2]; + bottom_right = u_test[((j + 1) / 2) * ctheta_int + (i - 1) / 2]; + } + else { + top_left = u_test[((j - 1) / 2) * ctheta_int]; + bottom_right = u_test[((j + 1) / 2) * ctheta_int + (i - 1) / 2]; + } + + double val = 0.5 * (top_left + bottom_right); + + EXPECT_NEAR(val, sol[j * p_level.ntheta_int + i], 1e-12) + << "Extrapolated Prolongation fails for Index (r,theta): (" + std::to_string(j) + "," + + std::to_string(i) + ")"; + ; + } + } + } +} + +INSTANTIATE_TEST_SUITE_P(Prolongation_size, test_prolongation, + ::testing::Values(std::make_tuple(0, false), std::make_tuple(1, false), + std::make_tuple(2, false), std::make_tuple(3, false), + std::make_tuple(4, false), std::make_tuple(5, false), + std::make_tuple(6, false), std::make_tuple(7, false), + std::make_tuple(8, false), std::make_tuple(0, true), + std::make_tuple(1, true), std::make_tuple(2, true), std::make_tuple(3, true), + std::make_tuple(4, true), std::make_tuple(5, true), std::make_tuple(6, true), + std::make_tuple(7, true), std::make_tuple(8, true))); \ No newline at end of file diff --git a/tests/test_restriction.cpp b/tests/test_restriction.cpp new file mode 100644 index 00000000..479d2dc8 --- /dev/null +++ b/tests/test_restriction.cpp @@ -0,0 +1,295 @@ +#include +#include "gmgpolar.h" +#include "mockgrid.h" +class test_restriction : public ::testing::TestWithParam> +{ +protected: + void SetUp() override + { + //initialize default parameters. + gyro::init_params(); + gyro::icntl[Param::verbose] = 0; + gyro::dcntl[Param::R0] = 1e-5; + gyro::f_grid_r = ""; + gyro::f_grid_theta = ""; + gyro::f_sol_in = ""; + gyro::f_sol_out = ""; + gyro::select_functions_class(gyro::icntl[Param::alpha_coeff], gyro::icntl[Param::beta_coeff], + gyro::icntl[Param::mod_pk], gyro::icntl[Param::prob]); + } +}; +/*! + * \brief Test the bilinear restriction operator used in the multigrid cycle coarse-grid correction. + * + * The test creates an arbitrary grid-function on the finer level and restricts it to the coarser level. + * On the coarse level we iterate over all nodes and accumulate the adjacent fine node values. + * This is matrix-free but it corresponds to applying the transposed prolongation matrix + * + * Parametrized tests are used to test for different grid sizes and with or without Dirichlet boundary conditions. + */ + +TEST_P(test_restriction, test_bilinear_restriction) +{ + gyro::icntl[Param::nr_exp] = (int)(std::get<0>(GetParam()) / 3) + 3; + gyro::icntl[Param::ntheta_exp] = (std::get<0>(GetParam()) % 3) + 3; + gyro::icntl[Param::DirBC_Interior] = std::get<1>(GetParam()); + + gmgpolar test_p; + + create_grid(test_p); //Mocking create_grid_polar, check_geom and define_coarse_nodes + + level& p_level = *(test_p.v_level[0]); + int ctheta_int = test_p.v_level[1]->ntheta_int; // number of coarse nodes in theta direction + int cr_int = test_p.v_level[1]->nr_int; + + p_level.m = test_p.v_level[0]->nr * test_p.v_level[0]->ntheta; //fine grid size + p_level.mc = test_p.v_level[1]->nr * test_p.v_level[1]->ntheta; //coarser grid size + + std::vector u_test(p_level.m); + for (int z = 0; z < p_level.m; z++) { + u_test[z] = + 1 - z + + pow(PI, + -z * z); //constructing arbitrary grid-function on coarse level to test our prolongation operator with. + } + + std::vector sol = p_level.apply_restriction_bi(u_test); + + /* + * The restriction operator accumulates the values in the direct vicinity and stores them in the corresponding node. + * These values are the eight fine nodes' value in the vicinity that are to be accumulated in the coarse node. + * + * We treat values in r as the x-axis. values in theta as the y-axis. Hence the vector 'adjacent' stores the values in the order: + * bottom_left, left, top_left, bottom, top, bottom_right, right, top_right, i.e. + * ^ --- 2 --- 4 --- 7 + * | --- 1 --- x --- 6 + * theta --- 0 --- 3 --- 5 + * | ------- r ------> + */ + + for (int j = 0; j < cr_int + 1; j++) { + for (int i = 0; i < ctheta_int; i++) { + std::vector adjacent(8, -1.0); + if (j == 0) { // innermost circle: nodes to the left (lower in radii indices) disappear + adjacent[0] = 0.0; + adjacent[1] = 0.0; + adjacent[2] = 0.0; + } + if (j == cr_int) { // outermost circle: nodes to the right disappear + adjacent[5] = 0.0; + adjacent[6] = 0.0; + adjacent[7] = 0.0; + } + + int k = 0; + // Store h_j,h_j-1 and k_i, k_i-1 for all adjacent nodes in the order as given above + std::vector> h_p(8, {-1, -1}); // (h_p, h_{p-1}) relative grid sizes + std::vector> k_q(8, {-1, -1}); // (k_q, k_{q-1}) + // z and y represent relative positions to the coarse node. e.g. if (z,y)=(-1,1) then we consider the fine node to the top-left + for (int z = -1; z < 2; z++) { + for (int y = -1; y < 2; y++) { + if (z != 0 || y != 0) { + if (adjacent[k] != 0.0) { + + adjacent[k] = + u_test[(2 * j + z) * p_level.ntheta_int + + ((i != 0 || y > -1) + ? (2 * i + y) + : p_level.ntheta_int - + 1)]; //adjacent value. consider periodic boundary in theta direction + + h_p[k] = {p_level.r[2 * j + z + 1] - p_level.r[2 * j + z], + p_level.r[2 * j + z] - + p_level.r[2 * j + z - 1]}; //distance in r coordinate for the adjacent node + + //to calculate k_q,k_{q-1} we consider the special case i=0 (2*i+y=y) separately. + if (i > 0) { + double indx = + (2 * i + y < p_level.ntheta_int - 1) ? p_level.theta[2 * i + y + 1] : 2 * PI; + + k_q[k] = {indx - p_level.theta[2 * i + y], + p_level.theta[2 * i + y] - p_level.theta[2 * i + y - 1]}; + } + else { + if (y == + 1) { //based on the value of y and periodic boundary conditions we get different values for {k_q,k_{q-1}} + k_q[k] = {p_level.theta[2] - p_level.theta[1], p_level.theta[1] - p_level.theta[0]}; + } + else if (y == 0) { + k_q[k] = {p_level.theta[1] - p_level.theta[0], + 2 * PI - p_level.theta[p_level.ntheta_int - 1]}; + } + else { + k_q[k] = {2 * PI - p_level.theta[p_level.ntheta_int - 1], + p_level.theta[p_level.ntheta_int - 1] - + p_level.theta[p_level.ntheta_int - 2]}; + } + } + } + k += 1; + } + } + } + /*values given by the operator. We multiply this with the grid-function value of the corresponding adjacent node*/ + + std::vector coeffs_hk{ + std::get<1>(h_p[0]) * std::get<1>(k_q[0]) / + ((std::get<0>(h_p[0]) + std::get<1>(h_p[0])) * (std::get<0>(k_q[0]) + std::get<1>(k_q[0]))), + std::get<1>(h_p[1]) / (std::get<0>(h_p[1]) + std::get<1>(h_p[1])), + std::get<1>(h_p[2]) * std::get<0>(k_q[2]) / + ((std::get<0>(h_p[2]) + std::get<1>(h_p[2])) * (std::get<0>(k_q[2]) + std::get<1>(k_q[2]))), + std::get<1>(k_q[3]) / (std::get<0>(k_q[3]) + std::get<1>(k_q[3])), + std::get<0>(k_q[4]) / (std::get<0>(k_q[4]) + std::get<1>(k_q[4])), + std::get<0>(h_p[5]) * std::get<1>(k_q[5]) / + ((std::get<0>(h_p[5]) + std::get<1>(h_p[5])) * (std::get<0>(k_q[5]) + std::get<1>(k_q[5]))), + std::get<0>(h_p[6]) / (std::get<0>(h_p[6]) + std::get<1>(h_p[6])), + std::get<0>(h_p[7]) * std::get<0>(k_q[7]) / + ((std::get<0>(h_p[7]) + std::get<1>(h_p[7])) * (std::get<0>(k_q[7]) + std::get<1>(k_q[7])))}; + + double finval = u_test[(2 * j) * p_level.ntheta_int + (2 * i)]; + for (int z = 0; z < 8; z++) { + finval += coeffs_hk[z] * adjacent[z]; //accumulate all values in the coarse node + } + EXPECT_NEAR(finval, sol[j * ctheta_int + i], 1e-10) + << "The test fails at Index for (r,theta): (" + std::to_string(j) + "," + std::to_string(i) + ")"; + } + } +} + +/*! + * \brief Test the injection restriction operator used in the implicit extrapolation step of the multigrid cycle. + * + * The Test creates an arbitrary grid-function on the finer level and restricts it. + * In this case this just corresponds to projecting the grid-function onto the coarse level. + * + * Parametrized tests are used to test for different grid sizes and with or without Dirichlet boundary conditions. + */ + +TEST_P(test_restriction, test_injection_restriction) +{ + gyro::icntl[Param::nr_exp] = (int)(std::get<0>(GetParam()) / 3) + 3; + gyro::icntl[Param::ntheta_exp] = (std::get<0>(GetParam()) % 3) + 3; + gyro::icntl[Param::DirBC_Interior] = std::get<1>(GetParam()); + + gmgpolar test_p; + + create_grid(test_p); + + level& p_level = *(test_p.v_level[0]); + int ctheta_int = test_p.v_level[1]->ntheta_int; + int cr_int = test_p.v_level[1]->nr_int; + + p_level.m = test_p.v_level[0]->nr * test_p.v_level[0]->ntheta; //number of nodes on fine level + p_level.mc = test_p.v_level[1]->nr * test_p.v_level[1]->ntheta; //number of nodes on coarse level + + std::vector u_test(p_level.m); + for (int z = 0; z < p_level.m; z++) { + u_test[z] = + 1 - z + + pow(PI, + -z * z); //constructing arbitrary grid-function on coarse level to test our prolongation operator with. + } + + std::vector sol = p_level.apply_restriction_inj(u_test); + + EXPECT_EQ((int)sol.size(), p_level.mc); + + for (int j = 0; j < cr_int + 1; j++) { + for (int i = 0; i < ctheta_int; i++) { + EXPECT_EQ(sol[j * ctheta_int + i], u_test[(2 * j) * p_level.ntheta_int + (2 * i)]) + << "The injection restriction fails at Index (r,theta): (" + std::to_string(j) + "," + + std::to_string(i) + ")"; //test all values + } + } +} + +/*! + * \brief Test the extrapolation restriction operator used in the implicit extrapolation step of the multigrid cycle. + * + * The Test creates an arbitrary grid-function on the finer level and restricts it. + * On the coarse level we iterate over all nodes and accumulate the adjacent fine node values based on the triangulation. + * We hence only consider 6 adjacent fine nodes, which lie on one of the edges spanned by the corresponding coarse node. + * This is matrix-free but it corresponds to applying the transposed prolongation matrix + * + * Parametrized tests are used to test for different grid sizes and with or without Dirichlet boundary conditions. + */ + +TEST_P(test_restriction, test_extrapolation_restriction) +{ + gyro::icntl[Param::nr_exp] = (int)(std::get<0>(GetParam()) / 3) + 3; + gyro::icntl[Param::ntheta_exp] = (std::get<0>(GetParam()) % 3) + 3; + gyro::icntl[Param::DirBC_Interior] = std::get<1>(GetParam()); + + gmgpolar test_p; + + create_grid(test_p); + + level& p_level = *(test_p.v_level[0]); + int ctheta_int = test_p.v_level[1]->ntheta_int; + int cr_int = test_p.v_level[1]->nr_int; + + p_level.m = test_p.v_level[0]->nr * test_p.v_level[0]->ntheta; //number of nodes on fine level + p_level.mc = test_p.v_level[1]->nr * test_p.v_level[1]->ntheta; //number of nodes on coarse level + + std::vector u_test(p_level.m); + for (int z = 0; z < p_level.mc; z++) { + u_test[z] = + 1 - z + + pow(PI, + -z * z); //constructing arbitrary grid-function on coarse level to test our prolongation operator with. + } + + std::vector sol = p_level.apply_restriction_ex(u_test); + + /* based on the triangulation, at most 6 adjacent fine nodes are considered: + * **left,top_left,bottom,top,bottom_right,right** + */ + for (int j = 0; j < cr_int + 1; j++) { + for (int i = 0; i < ctheta_int; i++) { + std::vector adjacent(6, -1.0); + if (j == 0) { // innermost circle + adjacent[0] = 0.0; + adjacent[1] = 0.0; + } + if (j == cr_int) { // outermost circle + adjacent[4] = 0.0; + adjacent[5] = 0.0; + } + + int k = 0; + //relative values of the fine nodes as in bilinear restriction + for (int z = -1; z < 2; z++) { + for (int y = -1; y < 2; y++) { + if ((z != 0 || y != 0) && (z * y != 1)) { + if (adjacent[k] != 0.0) { + + adjacent[k] = u_test[(2 * j + z) * p_level.ntheta_int + + ((i != 0 || y > -1) ? (2 * i + y) : p_level.ntheta_int - 1)]; + } + k += 1; + } + } + } + + double finval = u_test[(2 * j) * p_level.ntheta_int + (2 * i)]; + for (int z = 0; z < 6; z++) { + finval += + 0.5 * + adjacent[z]; // accumulate the values. the vector "coeffs_hk" reduces to 1/2 for every adjacent fine node + } + + EXPECT_NEAR(finval, sol[j * ctheta_int + i], 1e-12) + << "The test fails at Index for (r,theta): (" + std::to_string(j) + "," + std::to_string(i) + ")"; + } + } +} + +INSTANTIATE_TEST_SUITE_P(Restriction_size, test_restriction, + ::testing::Values(std::make_tuple(0, false), std::make_tuple(1, false), + std::make_tuple(2, false), std::make_tuple(3, false), + std::make_tuple(4, false), std::make_tuple(5, false), + std::make_tuple(6, false), std::make_tuple(7, false), + std::make_tuple(8, false), std::make_tuple(0, true), + std::make_tuple(1, true), std::make_tuple(2, true), std::make_tuple(3, true), + std::make_tuple(4, true), std::make_tuple(5, true), std::make_tuple(6, true), + std::make_tuple(7, true), std::make_tuple(8, true))); \ No newline at end of file