diff --git a/include/ExtrapolatedSmoother/ExtrapolatedSmootherGive/extrapolatedSmootherGive.h b/include/ExtrapolatedSmoother/ExtrapolatedSmootherGive/extrapolatedSmootherGive.h index 61aafb9d..981ba6dd 100644 --- a/include/ExtrapolatedSmoother/ExtrapolatedSmootherGive/extrapolatedSmootherGive.h +++ b/include/ExtrapolatedSmoother/ExtrapolatedSmootherGive/extrapolatedSmootherGive.h @@ -34,12 +34,14 @@ class ExtrapolatedSmootherGive : public ExtrapolatedSmoother // - 'radial_tridiagonal_solver_[index] refers to the radial line i_theta = 2*index, // - 'radial_diagonal_solver_[index] refers to the radial line i_theta = 2*index + 1. #ifdef GMGPOLAR_USE_MUMPS - SparseMatrixCOO inner_boundary_circle_matrix_; + using MatrixType = SparseMatrixCOO; DMUMPS_STRUC_C inner_boundary_mumps_solver_; #else - SparseMatrixCSR inner_boundary_circle_matrix_; + using MatrixType = SparseMatrixCSR; SparseLUSolver inner_boundary_lu_solver_; #endif + MatrixType inner_boundary_circle_matrix_; + std::vector> circle_diagonal_solver_; std::vector> radial_diagonal_solver_; std::vector> circle_tridiagonal_solver_; @@ -82,4 +84,12 @@ class ExtrapolatedSmootherGive : public ExtrapolatedSmoother void initializeMumpsSolver(DMUMPS_STRUC_C& mumps_solver, SparseMatrixCOO& solver_matrix); void finalizeMumpsSolver(DMUMPS_STRUC_C& mumps_solver); #endif + + void nodeBuildSmootherGive(int i_r, int i_theta, const PolarGrid& grid, bool DirBC_Interior, + MatrixType& inner_boundary_circle_matrix, + std::vector>& circle_diagonal_solver, + std::vector>& radial_diagonal_solver, + std::vector>& circle_tridiagonal_solver, + std::vector>& radial_tridiagonal_solver, double arr, + double att, double art, double detDF, double coeff_beta); }; diff --git a/src/ExtrapolatedSmoother/ExtrapolatedSmootherGive/buildAscMatrices.cpp b/src/ExtrapolatedSmoother/ExtrapolatedSmootherGive/buildAscMatrices.cpp index e6b72584..cabda889 100644 --- a/src/ExtrapolatedSmoother/ExtrapolatedSmootherGive/buildAscMatrices.cpp +++ b/src/ExtrapolatedSmoother/ExtrapolatedSmootherGive/buildAscMatrices.cpp @@ -3,1228 +3,1234 @@ #include "../../../include/common/geometry_helper.h" /* Tridiagonal matrices */ -#define UPDATE_TRIDIAGONAL_ELEMENT(matrix, row, column, value) \ - do { \ - if (row == column) \ - matrix.main_diagonal(row) += value; \ - else if (row == column - 1) \ - matrix.sub_diagonal(row) += value; \ - else if (row == 0 && column == matrix.columns() - 1) \ - matrix.cyclic_corner_element() += value; \ - } while (0) +inline void updateTridiagonalElement(SymmetricTridiagonalSolver& matrix, int row, int column, double value) +{ + if (row == column) + matrix.main_diagonal(row) += value; + else if (row == column - 1) + matrix.sub_diagonal(row) += value; + else if (row == 0 && column == matrix.columns() - 1) + matrix.cyclic_corner_element() += value; +} /* Diagonal matrices */ -#define UPDATE_DIAGONAL_ELEMENT(matrix, row, column, value) \ - do { \ - matrix.diagonal(row) += value; \ - } while (0) +inline void updateDiagonalElement(DiagonalSolver& matrix, int row, int column, double value) +{ + matrix.diagonal(row) += value; +} /* Inner Boundary COO/CSR matrix */ #ifdef GMGPOLAR_USE_MUMPS - #define COO_CSR_UPDATE(matrix, ptr, offset, row, col, val) \ - do { \ - matrix.row_index(ptr + offset) = row; \ - matrix.col_index(ptr + offset) = col; \ - matrix.value(ptr + offset) += val; \ - } while (0) +inline void updateCOOCSRMatrixElement(SparseMatrixCOO& matrix, int ptr, int offset, int row, int col, + double val) +{ + matrix.row_index(ptr + offset) = row; + matrix.col_index(ptr + offset) = col; + matrix.value(ptr + offset) += val; +} #else - #define COO_CSR_UPDATE(matrix, ptr, offset, row, col, val) \ - do { \ - matrix.row_nz_index(row, offset) = col; \ - matrix.row_nz_entry(row, offset) += val; \ - } while (0) +inline void updateCOOCSRMatrixElement(SparseMatrixCSR& matrix, int ptr, int offset, int row, int col, + double val) +{ + matrix.row_nz_index(row, offset) = col; + matrix.row_nz_entry(row, offset) += val; +} #endif -#define NODE_BUILD_SMOOTHER_GIVE(i_r, i_theta, grid, DirBC_Interior, inner_boundary_circle_matrix, \ - circle_diagonal_solver, circle_tridiagonal_solver, radial_diagonal_solver, \ - radial_tridiagonal_solver) \ - do { \ - assert(i_r >= 0 && i_r < grid.nr()); \ - assert(i_theta >= 0 && i_theta < grid.ntheta()); \ - \ - const int numberSmootherCircles = grid.numberSmootherCircles(); \ - const int lengthSmootherRadial = grid.lengthSmootherRadial(); \ - \ - assert(numberSmootherCircles >= 3); \ - assert(lengthSmootherRadial >= 3); \ - \ - int ptr, offset; \ - int row, column, col; \ - double value, val; \ - /* ------------------------------------------ */ \ - /* Node in the interior of the Circle Section */ \ - /* ------------------------------------------ */ \ - if (i_r > 0 && i_r < numberSmootherCircles - 1) { \ - double h1 = grid.radialSpacing(i_r - 1); \ - double h2 = grid.radialSpacing(i_r); \ - double k1 = grid.angularSpacing(i_theta - 1); \ - double k2 = grid.angularSpacing(i_theta); \ - double coeff1 = 0.5 * (k1 + k2) / h1; \ - double coeff2 = 0.5 * (k1 + k2) / h2; \ - double coeff3 = 0.5 * (h1 + h2) / k1; \ - double coeff4 = 0.5 * (h1 + h2) / k2; \ - \ - int i_theta_M1 = grid.wrapThetaIndex(i_theta - 1); \ - int i_theta_P1 = grid.wrapThetaIndex(i_theta + 1); \ - \ - int center_index = i_theta; \ - int left_index = i_theta; \ - int right_index = i_theta; \ - int bottom_index = i_theta_M1; \ - int top_index = i_theta_P1; \ - /* -------------------------- */ \ - /* Cyclic Tridiagonal Section */ \ - /* i_r % 2 == 1 */ \ - if (i_r & 1) { \ - /* i_theta % 2 == 1 */ \ - /* | x | o | x | */ \ - /* | | | | */ \ - /* | o | O | o | */ \ - /* | | | | */ \ - /* | x | o | x | */ \ - /* or */ \ - /* i_theta % 2 == 0 */ \ - /* | o | o | o | */ \ - /* | | | | */ \ - /* | x | O | x | */ \ - /* | | | | */ \ - /* | o | o | o | */ \ - \ - auto& center_matrix = circle_tridiagonal_solver[i_r / 2]; \ - auto& left_matrix = circle_diagonal_solver[(i_r - 1) / 2]; \ - auto& right_matrix = circle_diagonal_solver[(i_r + 1) / 2]; \ - \ - /* Fill matrix row of (i,j) */ \ - row = center_index; \ - column = center_index; \ - value = 0.25 * (h1 + h2) * (k1 + k2) * coeff_beta * fabs(detDF); /* Center: beta_{i,j} */ \ - UPDATE_TRIDIAGONAL_ELEMENT(center_matrix, row, column, value); \ - \ - row = center_index; \ - column = bottom_index; \ - value = -coeff3 * att; /* Bottom */ \ - UPDATE_TRIDIAGONAL_ELEMENT(center_matrix, row, column, value); \ - \ - row = center_index; \ - column = top_index; \ - value = -coeff4 * att; /* Top */ \ - UPDATE_TRIDIAGONAL_ELEMENT(center_matrix, row, column, value); \ - \ - row = center_index; \ - column = center_index; \ - value = (coeff1 + coeff2) * arr + (coeff3 + coeff4) * att; /* Center: (Left, Right, Bottom, Top) */ \ - UPDATE_TRIDIAGONAL_ELEMENT(center_matrix, row, column, value); \ - \ - /* Fill matrix row of (i,j-1) */ \ - row = bottom_index; \ - column = center_index; \ - value = -coeff3 * att; /* Top */ \ - UPDATE_TRIDIAGONAL_ELEMENT(center_matrix, row, column, value); \ - \ - row = bottom_index; \ - column = bottom_index; \ - value = coeff3 * att; /* Center: (Top) */ \ - UPDATE_TRIDIAGONAL_ELEMENT(center_matrix, row, column, value); \ - \ - /* Fill matrix row of (i,j+1) */ \ - row = top_index; \ - column = center_index; \ - value = -coeff4 * att; /* Bottom */ \ - UPDATE_TRIDIAGONAL_ELEMENT(center_matrix, row, column, value); \ - \ - row = top_index; \ - column = top_index; \ - value = coeff4 * att; /* Center: (Bottom) */ \ - UPDATE_TRIDIAGONAL_ELEMENT(center_matrix, row, column, value); \ - \ - if (i_theta & 1) { \ - /* i_theta % 2 == 1 */ \ - /* | x | o | x | */ \ - /* | | | | */ \ - /* | o | O | o | */ \ - /* | | | | */ \ - /* | x | o | x | */ \ - \ - /* Fill matrix row of (i-1,j) */ \ - if (i_r == 1) { \ - /* Only in the case of AcrossOrigin */ \ - if (!DirBC_Interior) { \ - row = left_index; \ - ptr = getCircleAscIndex(i_r - 1, i_theta); \ - \ - const Stencil& LeftStencil = getStencil(i_r - 1, i_theta); \ - \ - offset = LeftStencil[StencilPosition::Center]; \ - col = left_index; \ - val = +coeff1 * arr; /* Center: (Right) */ \ - COO_CSR_UPDATE(inner_boundary_circle_matrix, ptr, offset, row, col, val); \ - } \ - } \ - else { \ - row = left_index; \ - column = left_index; \ - value = coeff1 * arr; /* Center: (Right) */ \ - UPDATE_DIAGONAL_ELEMENT(left_matrix, row, column, value); \ - } \ - \ - /* Fill matrix row of (i+1,j) */ \ - row = right_index; \ - column = right_index; \ - value = coeff2 * arr; /* Center: (Left) */ \ - UPDATE_DIAGONAL_ELEMENT(right_matrix, row, column, value); \ - } \ - } \ - /* ---------------- */ \ - /* Diagonal Section */ \ - /* i_r % 2 == 0 */ \ - else { \ - /* i_theta % 2 == 1 */ \ - /* | o | x | o | */ \ - /* | | | | */ \ - /* | o | O | o | */ \ - /* | | | | */ \ - /* | o | x | o | */ \ - /* or */ \ - /* i_theta % 2 == 0 */ \ - /* | o | o | o | */ \ - /* | | | | */ \ - /* | o | X | o | */ \ - /* | | | | */ \ - /* | o | o | o | */ \ - \ - auto& center_matrix = circle_diagonal_solver[i_r / 2]; \ - auto& left_matrix = circle_tridiagonal_solver[(i_r - 1) / 2]; \ - auto& right_matrix = circle_tridiagonal_solver[(i_r + 1) / 2]; \ - \ - if (i_theta & 1) { /* i_theta % 2 == 1 */ \ - /* Fill matrix row of (i,j) */ \ - row = center_index; \ - column = center_index; \ - value = 0.25 * (h1 + h2) * (k1 + k2) * coeff_beta * fabs(detDF); /* Center: beta_{i,j} */ \ - UPDATE_DIAGONAL_ELEMENT(center_matrix, row, column, value); \ - \ - row = center_index; \ - column = center_index; \ - value = \ - (coeff1 + coeff2) * arr + (coeff3 + coeff4) * att; /* Center: (Left, Right, Bottom, Top) */ \ - UPDATE_DIAGONAL_ELEMENT(center_matrix, row, column, value); \ - } \ - else { /* i_theta % 2 == 0 */ \ - /* Fill matrix row of (i,j) */ \ - row = center_index; \ - column = center_index; \ - value = 1.0; \ - UPDATE_DIAGONAL_ELEMENT(center_matrix, row, column, value); \ - \ - /* Fill matrix row of (i,j-1) */ \ - row = bottom_index; \ - column = bottom_index; \ - value = coeff3 * att; /* Center: (Top) */ \ - UPDATE_DIAGONAL_ELEMENT(center_matrix, row, column, value); \ - \ - /* Fill matrix row of (i,j+1) */ \ - row = top_index; \ - column = top_index; \ - value = coeff4 * att; /* Center: (Bottom) */ \ - UPDATE_DIAGONAL_ELEMENT(center_matrix, row, column, value); \ - } \ - /* Fill matrix row of (i-1,j) */ \ - row = left_index; \ - column = left_index; \ - value = coeff1 * arr; /* Center: (Right) */ \ - UPDATE_TRIDIAGONAL_ELEMENT(left_matrix, row, column, value); \ - \ - /* Fill matrix row of (i+1,j) */ \ - row = right_index; \ - column = right_index; \ - value = coeff2 * arr; /* Center: (Left) */ \ - UPDATE_TRIDIAGONAL_ELEMENT(right_matrix, row, column, value); \ - } \ - } \ - /* ------------------------------------------ */ \ - /* Node in the interior of the Radial Section */ \ - /* ------------------------------------------ */ \ - else if (i_r > numberSmootherCircles && i_r < grid.nr() - 2) { \ - double h1 = grid.radialSpacing(i_r - 1); \ - double h2 = grid.radialSpacing(i_r); \ - double k1 = grid.angularSpacing(i_theta - 1); \ - double k2 = grid.angularSpacing(i_theta); \ - double coeff1 = 0.5 * (k1 + k2) / h1; \ - double coeff2 = 0.5 * (k1 + k2) / h2; \ - double coeff3 = 0.5 * (h1 + h2) / k1; \ - double coeff4 = 0.5 * (h1 + h2) / k2; \ - \ - int i_theta_M1 = grid.wrapThetaIndex(i_theta - 1); \ - int i_theta_P1 = grid.wrapThetaIndex(i_theta + 1); \ - \ - int center_index = i_r - numberSmootherCircles; \ - int left_index = i_r - numberSmootherCircles - 1; \ - int right_index = i_r - numberSmootherCircles + 1; \ - int bottom_index = i_r - numberSmootherCircles; \ - int top_index = i_r - numberSmootherCircles; \ - /* ------------------- */ \ - /* Tridiagonal Section */ \ - /* i_theta % 2 == 1 */ \ - if (i_theta & 1) { \ - /* i_r % 2 == 1 */ \ - /* ---------- */ \ - /* x o x */ \ - /* ---------- */ \ - /* o O o */ \ - /* ---------- */ \ - /* x o x */ \ - /* ---------- */ \ - /* or */ \ - /* i_r % 2 == 0 */ \ - /* ---------- */ \ - /* o x o */ \ - /* ---------- */ \ - /* o O o */ \ - /* ---------- */ \ - /* o x o */ \ - /* ---------- */ \ - \ - auto& center_matrix = radial_tridiagonal_solver[i_theta / 2]; \ - auto& bottom_matrix = radial_diagonal_solver[i_theta_M1 / 2]; \ - auto& top_matrix = radial_diagonal_solver[i_theta_P1 / 2]; \ - \ - /* Fill matrix row of (i,j) */ \ - row = center_index; \ - column = center_index; \ - value = 0.25 * (h1 + h2) * (k1 + k2) * coeff_beta * fabs(detDF); /* Center: beta_{i,j} */ \ - UPDATE_TRIDIAGONAL_ELEMENT(center_matrix, row, column, value); \ - \ - row = center_index; \ - column = left_index; \ - value = -coeff1 * arr; /* Left */ \ - UPDATE_TRIDIAGONAL_ELEMENT(center_matrix, row, column, value); \ - \ - row = center_index; \ - column = right_index; \ - value = -coeff2 * arr; /* Right */ \ - UPDATE_TRIDIAGONAL_ELEMENT(center_matrix, row, column, value); \ - \ - row = center_index; \ - column = center_index; \ - value = (coeff1 + coeff2) * arr + (coeff3 + coeff4) * att; /* Center: (Left, Right, Bottom, Top) */ \ - UPDATE_TRIDIAGONAL_ELEMENT(center_matrix, row, column, value); \ - \ - /* Fill matrix row of (i-1,j) */ \ - row = left_index; \ - column = center_index; \ - value = -coeff1 * arr; /* Right */ \ - UPDATE_TRIDIAGONAL_ELEMENT(center_matrix, row, column, value); \ - \ - row = left_index; \ - column = left_index; \ - value = coeff1 * arr; /* Center: (Right) */ \ - UPDATE_TRIDIAGONAL_ELEMENT(center_matrix, row, column, value); \ - \ - /* Fill matrix row of (i+1,j) */ \ - row = right_index; \ - column = center_index; \ - value = -coeff2 * arr; /* Left */ \ - UPDATE_TRIDIAGONAL_ELEMENT(center_matrix, row, column, value); \ - \ - row = right_index; \ - column = right_index; \ - value = coeff2 * arr; /* Center: (Left) */ \ - UPDATE_TRIDIAGONAL_ELEMENT(center_matrix, row, column, value); \ - \ - if (i_r & 1) { /* i_r % 2 == 1 */ \ - /* ---------- */ \ - /* x o x */ \ - /* ---------- */ \ - /* o O o */ \ - /* ---------- */ \ - /* x o x */ \ - /* ---------- */ \ - /* Fill matrix row of (i,j-1) */ \ - row = bottom_index; \ - column = bottom_index; \ - value = coeff3 * att; /* Center: (Top) */ \ - UPDATE_DIAGONAL_ELEMENT(bottom_matrix, row, column, value); \ - \ - /* Fill matrix row of (i,j+1) */ \ - row = top_index; \ - column = top_index; \ - value = coeff4 * att; /* Center: (Bottom) */ \ - UPDATE_DIAGONAL_ELEMENT(top_matrix, row, column, value); \ - } \ - } \ - /* ---------------- */ \ - /* Diagonal Section */ \ - /* i_theta % 2 == 0 */ \ - else { \ - /* i_r % 2 == 1 */ \ - /* ---------- */ \ - /* o o o */ \ - /* ---------- */ \ - /* x O x */ \ - /* ---------- */ \ - /* o o o */ \ - /* ---------- */ \ - /* or */ \ - /* i_r % 2 == 0 */ \ - /* ---------- */ \ - /* o o o */ \ - /* ---------- */ \ - /* o X o */ \ - /* ---------- */ \ - /* o o o */ \ - /* ---------- */ \ - \ - auto& center_matrix = radial_diagonal_solver[i_theta / 2]; \ - auto& bottom_matrix = radial_tridiagonal_solver[i_theta_M1 / 2]; \ - auto& top_matrix = radial_tridiagonal_solver[i_theta_P1 / 2]; \ - if (i_r & 1) { /* i_r % 2 == 1 */ \ - /* Fill matrix row of (i,j) */ \ - row = center_index; \ - column = center_index; \ - value = 0.25 * (h1 + h2) * (k1 + k2) * coeff_beta * fabs(detDF); /* Center: beta_{i,j} */ \ - UPDATE_DIAGONAL_ELEMENT(center_matrix, row, column, value); \ - \ - row = center_index; \ - column = center_index; \ - value = \ - (coeff1 + coeff2) * arr + (coeff3 + coeff4) * att; /* Center: (Left, Right, Bottom, Top) */ \ - UPDATE_DIAGONAL_ELEMENT(center_matrix, row, column, value); \ - } \ - else { /* i_r % 2 == 0 */ \ - /* Fill matrix row of (i,j) */ \ - row = center_index; \ - column = center_index; \ - value = 1.0; \ - UPDATE_DIAGONAL_ELEMENT(center_matrix, row, column, value); \ - \ - /* Fill matrix row of (i-1,j) */ \ - row = left_index; \ - column = left_index; \ - value = coeff1 * arr; /* Center: (Right) */ \ - UPDATE_DIAGONAL_ELEMENT(center_matrix, row, column, value); \ - \ - /* Fill matrix row of (i+1,j) */ \ - row = right_index; \ - column = right_index; \ - value = coeff2 * arr; /* Center: (Left) */ \ - UPDATE_DIAGONAL_ELEMENT(center_matrix, row, column, value); \ - } \ - /* Fill matrix row of (i,j-1) */ \ - row = bottom_index; \ - column = bottom_index; \ - value = coeff3 * att; /* Center: (Top) */ \ - UPDATE_TRIDIAGONAL_ELEMENT(bottom_matrix, row, column, value); \ - \ - /* Fill matrix row of (i,j+1) */ \ - row = top_index; \ - column = top_index; \ - value = coeff4 * att; /* Center: (Bottom) */ \ - UPDATE_TRIDIAGONAL_ELEMENT(top_matrix, row, column, value); \ - } \ - } \ - /* ------------------------------------------ */ \ - /* Circle Section: Node in the inner boundary */ \ - /* ------------------------------------------ */ \ - else if (i_r == 0) { \ - /* ------------------------------------------------ */ \ - /* Case 1: Dirichlet boundary on the inner boundary */ \ - /* ------------------------------------------------ */ \ - if (DirBC_Interior) { \ - /* Fill result(i,j) */ \ - double h2 = grid.radialSpacing(i_r); \ - double k1 = grid.angularSpacing(i_theta - 1); \ - double k2 = grid.angularSpacing(i_theta); \ - double coeff2 = 0.5 * (k1 + k2) / h2; \ - \ - int i_theta_M1 = grid.wrapThetaIndex(i_theta - 1); \ - int i_theta_P1 = grid.wrapThetaIndex(i_theta + 1); \ - \ - auto& center_matrix = inner_boundary_circle_matrix; \ - auto& right_matrix = circle_tridiagonal_solver[(i_r + 1) / 2]; \ - \ - int center_index = i_theta; \ - int right_index = i_theta; \ - int bottom_index = i_theta_M1; \ - int top_index = i_theta_P1; \ - \ - /* Fill matrix row of (i,j) */ \ - row = center_index; \ - ptr = getCircleAscIndex(i_r, i_theta); \ - \ - const Stencil& CenterStencil = getStencil(i_r, i_theta); \ - \ - offset = CenterStencil[StencilPosition::Center]; \ - col = center_index; \ - val = 1.0; \ - COO_CSR_UPDATE(center_matrix, ptr, offset, row, col, val); \ - \ - /* Fill matrix row of (i+1,j) */ \ - row = right_index; \ - column = right_index; \ - value = coeff2 * arr; /* Center: (Left) */ \ - UPDATE_TRIDIAGONAL_ELEMENT(right_matrix, row, column, value); \ - } \ - else { \ - /* ------------------------------------------------------------- */ \ - /* Case 2: Across origin discretization on the interior boundary */ \ - /* ------------------------------------------------------------- */ \ - /* h1 gets replaced with 2 * R0. */ \ - /* (i_r-1,i_theta) gets replaced with (i_r, i_theta + (grid.ntheta()/2)). */ \ - /* Some more adjustments from the changing the 9-point stencil to the artifical 7-point stencil. */ \ - double h1 = 2.0 * grid.radius(0); \ - double h2 = grid.radialSpacing(i_r); \ - double k1 = grid.angularSpacing(i_theta - 1); \ - double k2 = grid.angularSpacing(i_theta); \ - double coeff1 = 0.5 * (k1 + k2) / h1; \ - double coeff2 = 0.5 * (k1 + k2) / h2; \ - double coeff3 = 0.5 * (h1 + h2) / k1; \ - double coeff4 = 0.5 * (h1 + h2) / k2; \ - \ - const int i_theta_M1 = grid.wrapThetaIndex(i_theta - 1); \ - const int i_theta_P1 = grid.wrapThetaIndex(i_theta + 1); \ - const int i_theta_AcrossOrigin = grid.wrapThetaIndex(i_theta + (grid.ntheta() / 2)); \ - \ - const int center_index = i_theta; \ - const int left_index = i_theta_AcrossOrigin; \ - const int right_index = i_theta; \ - const int bottom_index = i_theta_M1; \ - const int top_index = i_theta_P1; \ - \ - const int center_nz_index = getCircleAscIndex(i_r, i_theta); \ - const int bottom_nz_index = getCircleAscIndex(i_r, i_theta_M1); \ - const int top_nz_index = getCircleAscIndex(i_r, i_theta_P1); \ - const int left_nz_index = getCircleAscIndex(i_r, i_theta_AcrossOrigin); \ - \ - int nz_index; \ - const Stencil& CenterStencil = getStencil(i_r, i_theta); \ - \ - if (i_theta & 1) { \ - /* i_theta % 2 == 1 */ \ - /* -| x | o | x | */ \ - /* -| | | | */ \ - /* -| O | o | o | */ \ - /* -| | | | */ \ - /* -| x | o | x | */ \ - \ - auto& center_matrix = inner_boundary_circle_matrix; \ - auto& right_matrix = circle_tridiagonal_solver[(i_r + 1) / 2]; \ - auto& left_matrix = inner_boundary_circle_matrix; \ - /* Fill matrix row of (i,j) */ \ - row = center_index; \ - ptr = center_nz_index; \ - \ - offset = CenterStencil[StencilPosition::Center]; \ - col = center_index; \ - val = 0.25 * (h1 + h2) * (k1 + k2) * coeff_beta * fabs(detDF); /* beta_{i,j} */ \ - COO_CSR_UPDATE(center_matrix, ptr, offset, row, col, val); \ - \ - offset = CenterStencil[StencilPosition::Left]; \ - col = left_index; \ - val = -coeff1 * arr; /* Left */ \ - COO_CSR_UPDATE(center_matrix, ptr, offset, row, col, val); \ - \ - offset = CenterStencil[StencilPosition::Center]; \ - col = center_index; \ - val = (coeff1 + coeff2) * arr + (coeff3 + coeff4) * att; /* Center: (Left, Right, Bottom, Top) */ \ - COO_CSR_UPDATE(center_matrix, ptr, offset, row, col, val); \ - \ - /* Fill matrix row of (i-1,j) */ \ - /* From view the view of the across origin node, */ \ - /* the directions are roatated by 180 degrees in the stencil! */ \ - row = left_index; \ - ptr = left_nz_index; \ - \ - const Stencil& LeftStencil = CenterStencil; \ - \ - offset = LeftStencil[StencilPosition::Left]; \ - col = center_index; \ - val = -coeff1 * arr; /* Right -> Left*/ \ - COO_CSR_UPDATE(left_matrix, ptr, offset, row, col, val); \ - \ - offset = LeftStencil[StencilPosition::Center]; \ - col = left_index; \ - val = +coeff1 * arr; /* Center: (Right) -> Center: (Left) */ \ - COO_CSR_UPDATE(left_matrix, ptr, offset, row, col, val); \ - \ - /* Fill matrix row of (i+1,j) */ \ - row = right_index; \ - column = right_index; \ - value = coeff2 * arr; /* Center: (Left) */ \ - UPDATE_TRIDIAGONAL_ELEMENT(right_matrix, row, column, value); \ - } \ - else { \ - /* i_theta % 2 == 0 */ \ - /* -| o | o | o | */ \ - /* -| | | | */ \ - /* -| X | o | x | */ \ - /* -| | | | */ \ - /* -| o | o | o | */ \ - \ - auto& center_matrix = inner_boundary_circle_matrix; \ - auto& right_matrix = circle_tridiagonal_solver[(i_r + 1) / 2]; \ - auto& left_matrix = inner_boundary_circle_matrix; \ - /* Fill matrix row of (i,j) */ \ - row = center_index; \ - ptr = center_nz_index; \ - \ - offset = CenterStencil[StencilPosition::Center]; \ - col = center_index; \ - val = 1.0; \ - COO_CSR_UPDATE(center_matrix, ptr, offset, row, col, val); \ - \ - /* Fill matrix row of (i,j-1) */ \ - row = bottom_index; \ - ptr = bottom_nz_index; \ - \ - const Stencil& BottomStencil = CenterStencil; \ - \ - offset = BottomStencil[StencilPosition::Center]; \ - col = bottom_index; \ - val = +coeff3 * att; /* Center: (Top) */ \ - COO_CSR_UPDATE(center_matrix, ptr, offset, row, col, val); \ - \ - /* Fill matrix row of (i,j+1) */ \ - row = top_index; \ - ptr = top_nz_index; \ - \ - const Stencil& TopStencil = CenterStencil; \ - \ - offset = TopStencil[StencilPosition::Center]; \ - col = top_index; \ - val = +coeff4 * att; /* Center: (Bottom) */ \ - COO_CSR_UPDATE(center_matrix, ptr, offset, row, col, val); \ - \ - /* Fill matrix row of (i+1,j) */ \ - row = right_index; \ - column = right_index; \ - value = coeff2 * arr; /* Center: (Left) */ \ - UPDATE_TRIDIAGONAL_ELEMENT(right_matrix, row, column, value); \ - } \ - } \ - } \ - /* ------------------------------------------- */ \ - /* Circle Section: Node next to radial section */ \ - /* ------------------------------------------- */ \ - else if (i_r == numberSmootherCircles - 1) { \ - assert(i_r > 1); \ - \ - double h1 = grid.radialSpacing(i_r - 1); \ - double h2 = grid.radialSpacing(i_r); \ - double k1 = grid.angularSpacing(i_theta - 1); \ - double k2 = grid.angularSpacing(i_theta); \ - double coeff1 = 0.5 * (k1 + k2) / h1; \ - double coeff2 = 0.5 * (k1 + k2) / h2; \ - double coeff3 = 0.5 * (h1 + h2) / k1; \ - double coeff4 = 0.5 * (h1 + h2) / k2; \ - \ - int i_theta_M1 = grid.wrapThetaIndex(i_theta - 1); \ - int i_theta_P1 = grid.wrapThetaIndex(i_theta + 1); \ - \ - int center_index = i_theta; \ - int left_index = i_theta; \ - int right_index = 0; \ - int bottom_index = i_theta_M1; \ - int top_index = i_theta_P1; \ - \ - if (i_r & 1) { \ - if (i_theta & 1) { \ - /* i_r % 2 == 1 and i_theta % 2 == 1 */ \ - /* | o | x | o || x o x o */ \ - /* | | | || -------------- */ \ - /* | o | o | O || o o o o */ \ - /* | | | || -------------- */ \ - /* | o | x | o || x o x o */ \ - \ - auto& center_matrix = circle_tridiagonal_solver[i_r / 2]; \ - auto& left_matrix = circle_diagonal_solver[(i_r - 1) / 2]; \ - auto& right_matrix = radial_tridiagonal_solver[i_theta / 2]; \ - \ - /* Fill matrix row of (i,j) */ \ - row = center_index; \ - column = center_index; \ - value = 0.25 * (h1 + h2) * (k1 + k2) * coeff_beta * fabs(detDF); /* Center: beta_{i,j} */ \ - UPDATE_TRIDIAGONAL_ELEMENT(center_matrix, row, column, value); \ - \ - row = center_index; \ - column = bottom_index; \ - value = -coeff3 * att; /* Bottom */ \ - UPDATE_TRIDIAGONAL_ELEMENT(center_matrix, row, column, value); \ - \ - row = center_index; \ - column = top_index; \ - value = -coeff4 * att; /* Top */ \ - UPDATE_TRIDIAGONAL_ELEMENT(center_matrix, row, column, value); \ - \ - row = center_index; \ - column = center_index; \ - value = \ - (coeff1 + coeff2) * arr + (coeff3 + coeff4) * att; /* Center: (Left, Right, Bottom, Top) */ \ - UPDATE_TRIDIAGONAL_ELEMENT(center_matrix, row, column, value); \ - \ - /* Fill matrix row of (i,j-1) */ \ - row = bottom_index; \ - column = center_index; \ - value = -coeff3 * att; /* Top */ \ - UPDATE_TRIDIAGONAL_ELEMENT(center_matrix, row, column, value); \ - \ - row = bottom_index; \ - column = bottom_index; \ - value = coeff3 * att; /* Center: (Top) */ \ - UPDATE_TRIDIAGONAL_ELEMENT(center_matrix, row, column, value); \ - \ - /* Fill matrix row of (i,j+1) */ \ - row = top_index; \ - column = center_index; \ - value = -coeff4 * att; /* Bottom */ \ - UPDATE_TRIDIAGONAL_ELEMENT(center_matrix, row, column, value); \ - \ - row = top_index; \ - column = top_index; \ - value = coeff4 * att; /* Center: (Bottom) */ \ - UPDATE_TRIDIAGONAL_ELEMENT(center_matrix, row, column, value); \ - \ - /* Fill matrix row of (i-1,j) */ \ - row = left_index; \ - column = left_index; \ - value = coeff1 * arr; /* Center: (Right) */ \ - UPDATE_DIAGONAL_ELEMENT(left_matrix, row, column, value); \ - \ - /* Fill matrix row of (i+1,j) */ \ - row = right_index; \ - column = right_index; \ - value = coeff2 * arr; /* Center: (Left) */ \ - UPDATE_TRIDIAGONAL_ELEMENT(right_matrix, row, column, value); \ - } \ - else { \ - /* i_r % 2 == 1 and i_theta % 2 == 0 */ \ - /* | o | o | o || o o o o */ \ - /* | | | || -------------- */ \ - /* | o | x | O || x o x o */ \ - /* | | | || -------------- */ \ - /* | o | o | o || o o o o */ \ - \ - auto& center_matrix = circle_tridiagonal_solver[i_r / 2]; \ - auto& left_matrix = circle_diagonal_solver[(i_r - 1) / 2]; \ - auto& right_matrix = radial_diagonal_solver[i_theta / 2]; \ - \ - /* Fill matrix row of (i,j) */ \ - row = center_index; \ - column = center_index; \ - value = 0.25 * (h1 + h2) * (k1 + k2) * coeff_beta * fabs(detDF); /* Center: beta_{i,j} */ \ - UPDATE_TRIDIAGONAL_ELEMENT(center_matrix, row, column, value); \ - \ - row = center_index; \ - column = bottom_index; \ - value = -coeff3 * att; /* Bottom */ \ - UPDATE_TRIDIAGONAL_ELEMENT(center_matrix, row, column, value); \ - \ - row = center_index; \ - column = top_index; \ - value = -coeff4 * att; /* Top */ \ - UPDATE_TRIDIAGONAL_ELEMENT(center_matrix, row, column, value); \ - \ - row = center_index; \ - column = center_index; \ - value = \ - (coeff1 + coeff2) * arr + (coeff3 + coeff4) * att; /* Center: (Left, Right, Bottom, Top) */ \ - UPDATE_TRIDIAGONAL_ELEMENT(center_matrix, row, column, value); \ - \ - /* Fill matrix row of (i,j-1) */ \ - row = bottom_index; \ - column = center_index; \ - value = -coeff3 * att; /* Top */ \ - UPDATE_TRIDIAGONAL_ELEMENT(center_matrix, row, column, value); \ - \ - row = bottom_index; \ - column = bottom_index; \ - value = coeff3 * att; /* Center: (Top) */ \ - UPDATE_TRIDIAGONAL_ELEMENT(center_matrix, row, column, value); \ - \ - /* Fill matrix row of (i,j+1) */ \ - row = top_index; \ - column = center_index; \ - value = -coeff4 * att; /* Bottom */ \ - UPDATE_TRIDIAGONAL_ELEMENT(center_matrix, row, column, value); \ - \ - row = top_index; \ - column = top_index; \ - value = coeff4 * att; /* Center: (Bottom) */ \ - UPDATE_TRIDIAGONAL_ELEMENT(center_matrix, row, column, value); \ - } \ - } \ - else { \ - if (i_theta & 1) { \ - /* i_r % 2 == 0 and i_theta % 2 == 1 */ \ - /* | x | o | x || o x o x */ \ - /* | | | || -------------- */ \ - /* | o | o | O || o o o o */ \ - /* | | | || -------------- */ \ - /* | x | o | x || o x o x */ \ - \ - auto& center_matrix = circle_diagonal_solver[i_r / 2]; \ - auto& left_matrix = circle_tridiagonal_solver[(i_r - 1) / 2]; \ - auto& right_matrix = radial_tridiagonal_solver[i_theta / 2]; \ - \ - /* Fill matrix row of (i,j) */ \ - row = center_index; \ - column = center_index; \ - value = 0.25 * (h1 + h2) * (k1 + k2) * coeff_beta * fabs(detDF); /* Center: beta_{i,j} */ \ - UPDATE_DIAGONAL_ELEMENT(center_matrix, row, column, value); \ - \ - row = center_index; \ - column = center_index; \ - value = \ - (coeff1 + coeff2) * arr + (coeff3 + coeff4) * att; /* Center: (Left, Right, Bottom, Top) */ \ - UPDATE_DIAGONAL_ELEMENT(center_matrix, row, column, value); \ - \ - /* Fill matrix row of (i-1,j) */ \ - row = left_index; \ - column = left_index; \ - value = coeff1 * arr; /* Center: (Right) */ \ - UPDATE_TRIDIAGONAL_ELEMENT(left_matrix, row, column, value); \ - \ - /* Fill matrix row of (i+1,j) */ \ - row = right_index; \ - column = right_index; \ - value = coeff2 * arr; /* Center: (Left) */ \ - UPDATE_TRIDIAGONAL_ELEMENT(right_matrix, row, column, value); \ - } \ - else { \ - /* i_r % 2 == 0 and i_theta % 2 == 0 */ \ - /* | o | o | o || o o o o */ \ - /* | | | || -------------- */ \ - /* | x | o | X || o x o x */ \ - /* | | | || -------------- */ \ - /* | o | o | o || o o o o */ \ - \ - auto& center_matrix = circle_diagonal_solver[i_r / 2]; \ - auto& left_matrix = circle_tridiagonal_solver[(i_r - 1) / 2]; \ - auto& right_matrix = radial_diagonal_solver[i_theta / 2]; \ - \ - /* Fill matrix row of (i,j) */ \ - row = center_index; \ - column = center_index; \ - value = 1.0; \ - UPDATE_DIAGONAL_ELEMENT(center_matrix, row, column, value); \ - \ - /* Fill matrix row of (i,j-1) */ \ - row = bottom_index; \ - column = bottom_index; \ - value = coeff3 * att; /* Center: (Top) */ \ - UPDATE_DIAGONAL_ELEMENT(center_matrix, row, column, value); \ - \ - /* Fill matrix row of (i,j+1) */ \ - row = top_index; \ - column = top_index; \ - value = coeff4 * att; /* Center: (Bottom) */ \ - UPDATE_DIAGONAL_ELEMENT(center_matrix, row, column, value); \ - \ - /* Fill matrix row of (i-1,j) */ \ - row = left_index; \ - column = left_index; \ - value = coeff1 * arr; /* Center: (Right) */ \ - UPDATE_TRIDIAGONAL_ELEMENT(left_matrix, row, column, value); \ - \ - /* Fill matrix row of (i+1,j) */ \ - row = right_index; \ - column = right_index; \ - value = coeff2 * arr; /* Center: (Left) */ \ - UPDATE_DIAGONAL_ELEMENT(right_matrix, row, column, value); \ - } \ - } \ - } \ - /* --------------------------------------------- */ \ - /* Radial Section: Node next to circular section */ \ - /* --------------------------------------------- */ \ - else if (i_r == numberSmootherCircles) { \ - double h1 = grid.radialSpacing(i_r - 1); \ - double h2 = grid.radialSpacing(i_r); \ - double k1 = grid.angularSpacing(i_theta - 1); \ - double k2 = grid.angularSpacing(i_theta); \ - double coeff1 = 0.5 * (k1 + k2) / h1; \ - double coeff2 = 0.5 * (k1 + k2) / h2; \ - double coeff3 = 0.5 * (h1 + h2) / k1; \ - double coeff4 = 0.5 * (h1 + h2) / k2; \ - \ - const int i_theta_M1 = grid.wrapThetaIndex(i_theta - 1); \ - const int i_theta_P1 = grid.wrapThetaIndex(i_theta + 1); \ - \ - const int center_index = i_r - numberSmootherCircles; \ - const int left_index = i_theta; \ - const int right_index = i_r - numberSmootherCircles + 1; \ - const int bottom_index = i_r - numberSmootherCircles; \ - const int top_index = i_r - numberSmootherCircles; \ - \ - if (i_theta & 1) { \ - if (i_r & 1) { \ - /* i_theta % 2 == 1 and i_r % 2 == 1 */ \ - /* | x | o | x || o x o x */ \ - /* | | | || -------------- */ \ - /* | o | o | o || O o o o */ \ - /* | | | || -------------- */ \ - /* | x | o | x || o x o x */ \ - \ - auto& center_matrix = radial_tridiagonal_solver[i_theta / 2]; \ - auto& bottom_matrix = radial_diagonal_solver[i_theta_M1 / 2]; \ - auto& top_matrix = radial_diagonal_solver[i_theta_P1 / 2]; \ - auto& left_matrix = circle_diagonal_solver[(i_r - 1) / 2]; \ - \ - /* Fill matrix row of (i,j) */ \ - row = center_index; \ - column = center_index; \ - value = 0.25 * (h1 + h2) * (k1 + k2) * coeff_beta * fabs(detDF); /* Center: beta_{i,j} */ \ - UPDATE_TRIDIAGONAL_ELEMENT(center_matrix, row, column, value); \ - \ - row = center_index; \ - column = right_index; \ - value = -coeff2 * arr; /* Right */ \ - UPDATE_TRIDIAGONAL_ELEMENT(center_matrix, row, column, value); \ - \ - row = center_index; \ - column = center_index; \ - value = \ - (coeff1 + coeff2) * arr + (coeff3 + coeff4) * att; /* Center: (Left, Right, Bottom, Top) */ \ - UPDATE_TRIDIAGONAL_ELEMENT(center_matrix, row, column, value); \ - \ - /* Fill matrix row of (i-1,j) */ \ - row = left_index; \ - column = left_index; \ - value = coeff1 * arr; /* Center: (Right) */ \ - UPDATE_DIAGONAL_ELEMENT(left_matrix, row, column, value); \ - \ - /* Fill matrix row of (i+1,j) */ \ - row = right_index; \ - column = center_index; \ - value = -coeff2 * arr; /* Left */ \ - UPDATE_TRIDIAGONAL_ELEMENT(center_matrix, row, column, value); \ - \ - row = right_index; \ - column = right_index; \ - value = coeff2 * arr; /* Center: (Left) */ \ - UPDATE_TRIDIAGONAL_ELEMENT(center_matrix, row, column, value); \ - \ - /* Fill matrix row of (i,j-1) */ \ - row = bottom_index; \ - column = bottom_index; \ - value = coeff3 * att; /* Center: (Top) */ \ - UPDATE_DIAGONAL_ELEMENT(bottom_matrix, row, column, value); \ - \ - /* Fill matrix row of (i,j+1) */ \ - row = top_index; \ - column = top_index; \ - value = coeff4 * att; /* Center: (Bottom) */ \ - UPDATE_DIAGONAL_ELEMENT(top_matrix, row, column, value); \ - } \ - else { \ - /* i_theta % 2 == 1 and i_r % 2 == 0 */ \ - /* | o | x | o || x o x o */ \ - /* | | | || -------------- */ \ - /* | o | o | o || O o o o */ \ - /* | | | || -------------- */ \ - /* | o | x | o || x o x o */ \ - \ - auto& center_matrix = radial_tridiagonal_solver[i_theta / 2]; \ - auto& bottom_matrix = radial_diagonal_solver[i_theta_M1 / 2]; \ - auto& top_matrix = radial_diagonal_solver[i_theta_P1 / 2]; \ - auto& left_matrix = circle_tridiagonal_solver[(i_r - 1) / 2]; \ - \ - /* Fill matrix row of (i,j) */ \ - row = center_index; \ - column = center_index; \ - value = 0.25 * (h1 + h2) * (k1 + k2) * coeff_beta * fabs(detDF); /* Center: beta_{i,j} */ \ - UPDATE_TRIDIAGONAL_ELEMENT(center_matrix, row, column, value); \ - \ - row = center_index; \ - column = right_index; \ - value = -coeff2 * arr; /* Right */ \ - UPDATE_TRIDIAGONAL_ELEMENT(center_matrix, row, column, value); \ - \ - row = center_index; \ - column = center_index; \ - value = \ - (coeff1 + coeff2) * arr + (coeff3 + coeff4) * att; /* Center: (Left, Right, Bottom, Top) */ \ - UPDATE_TRIDIAGONAL_ELEMENT(center_matrix, row, column, value); \ - \ - /* Fill matrix row of (i-1,j) */ \ - row = left_index; \ - column = left_index; \ - value = coeff1 * arr; /* Center: (Right) */ \ - UPDATE_TRIDIAGONAL_ELEMENT(left_matrix, row, column, value); \ - \ - /* Fill matrix row of (i+1,j) */ \ - row = right_index; \ - column = center_index; \ - value = -coeff2 * arr; /* Left */ \ - UPDATE_TRIDIAGONAL_ELEMENT(center_matrix, row, column, value); \ - \ - row = right_index; \ - column = right_index; \ - value = coeff2 * arr; /* Center: (Left) */ \ - UPDATE_TRIDIAGONAL_ELEMENT(center_matrix, row, column, value); \ - } \ - } \ - else { \ - if (i_r & 1) { \ - /* i_theta % 2 == 0 and i_r % 2 == 1 */ \ - /* | o | o | o || o o o o */ \ - /* | | | || -------------- */ \ - /* | x | o | x || O x o x */ \ - /* | | | || -------------- */ \ - /* | o | o | o || o o o o */ \ - \ - auto& center_matrix = radial_diagonal_solver[i_theta / 2]; \ - auto& bottom_matrix = radial_tridiagonal_solver[i_theta_M1 / 2]; \ - auto& top_matrix = radial_tridiagonal_solver[i_theta_P1 / 2]; \ - auto& left_matrix = circle_diagonal_solver[(i_r - 1) / 2]; \ - \ - /* Fill matrix row of (i,j) */ \ - row = center_index; \ - column = center_index; \ - value = 0.25 * (h1 + h2) * (k1 + k2) * coeff_beta * fabs(detDF); /* Center: beta_{i,j} */ \ - UPDATE_DIAGONAL_ELEMENT(center_matrix, row, column, value); \ - \ - row = center_index; \ - column = center_index; \ - value = \ - (coeff1 + coeff2) * arr + (coeff3 + coeff4) * att; /* Center: (Left, Right, Bottom, Top) */ \ - UPDATE_DIAGONAL_ELEMENT(center_matrix, row, column, value); \ - \ - /* Fill matrix row of (i,j-1) */ \ - row = bottom_index; \ - column = bottom_index; \ - value = coeff3 * att; /* Center: (Top) */ \ - UPDATE_TRIDIAGONAL_ELEMENT(bottom_matrix, row, column, value); \ - \ - /* Fill matrix row of (i,j+1) */ \ - row = top_index; \ - column = top_index; \ - value = coeff4 * att; /* Center: (Bottom) */ \ - UPDATE_TRIDIAGONAL_ELEMENT(top_matrix, row, column, value); \ - } \ - else { \ - /* i_theta % 2 == 0 and i_r % 2 == 0 */ \ - /* | o | o | o || o o o o */ \ - /* | | | || -------------- */ \ - /* | o | x | o || X o x o */ \ - /* | | | || -------------- */ \ - /* | o | o | o || o o o o */ \ - \ - auto& center_matrix = radial_diagonal_solver[i_theta / 2]; \ - auto& bottom_matrix = radial_tridiagonal_solver[i_theta_M1 / 2]; \ - auto& top_matrix = radial_tridiagonal_solver[i_theta_P1 / 2]; \ - auto& left_matrix = circle_tridiagonal_solver[(i_r - 1) / 2]; \ - \ - /* Fill matrix row of (i,j) */ \ - \ - row = center_index; \ - column = center_index; \ - value = 1.0; \ - UPDATE_DIAGONAL_ELEMENT(center_matrix, row, column, value); \ - \ - /* Fill matrix row of (i-1,j) */ \ - row = left_index; \ - column = left_index; \ - value = coeff1 * arr; /* Center: (Right) */ \ - UPDATE_TRIDIAGONAL_ELEMENT(left_matrix, row, column, value); \ - \ - /* Fill matrix row of (i+1,j) */ \ - row = right_index; \ - column = right_index; \ - value = coeff2 * arr; /* Center: (Left) */ \ - UPDATE_DIAGONAL_ELEMENT(center_matrix, row, column, value); \ - \ - /* Fill matrix row of (i,j-1) */ \ - row = bottom_index; \ - column = bottom_index; \ - value = coeff3 * att; /* Center: (Top) */ \ - UPDATE_TRIDIAGONAL_ELEMENT(bottom_matrix, row, column, value); \ - \ - /* Fill matrix row of (i,j+1) */ \ - row = top_index; \ - column = top_index; \ - value = coeff4 * att; /* Center: (Bottom) */ \ - UPDATE_TRIDIAGONAL_ELEMENT(top_matrix, row, column, value); \ - } \ - } \ - } \ - /* ------------------------------------------- */ \ - /* Radial Section: Node next to outer boundary */ \ - /* ------------------------------------------- */ \ - else if (i_r == grid.nr() - 2) { \ - assert(i_r % 2 == 1); \ - \ - double h1 = grid.radialSpacing(i_r - 1); \ - double h2 = grid.radialSpacing(i_r); \ - double k1 = grid.angularSpacing(i_theta - 1); \ - double k2 = grid.angularSpacing(i_theta); \ - double coeff1 = 0.5 * (k1 + k2) / h1; \ - double coeff2 = 0.5 * (k1 + k2) / h2; \ - double coeff3 = 0.5 * (h1 + h2) / k1; \ - double coeff4 = 0.5 * (h1 + h2) / k2; \ - \ - int i_theta_M1 = grid.wrapThetaIndex(i_theta - 1); \ - int i_theta_P1 = grid.wrapThetaIndex(i_theta + 1); \ - \ - int center_index = i_r - numberSmootherCircles; \ - int left_index = i_r - numberSmootherCircles - 1; \ - int right_index = i_r - numberSmootherCircles + 1; \ - int bottom_index = i_r - numberSmootherCircles; \ - int top_index = i_r - numberSmootherCircles; \ - \ - if (i_theta & 1) { \ - /* i_theta % 2 == 1 */ \ - /* ---------------|| */ \ - /* o x o x || */ \ - /* ---------------|| */ \ - /* o o O o || */ \ - /* ---------------|| */ \ - /* o x o x || */ \ - /* ---------------|| */ \ - auto& center_matrix = radial_tridiagonal_solver[i_theta / 2]; \ - auto& bottom_matrix = radial_diagonal_solver[i_theta_M1 / 2]; \ - auto& top_matrix = radial_diagonal_solver[i_theta_P1 / 2]; \ - /* Fill matrix row of (i,j) */ \ - row = center_index; \ - column = center_index; \ - value = 0.25 * (h1 + h2) * (k1 + k2) * coeff_beta * fabs(detDF); /* Center: beta_{i,j} */ \ - UPDATE_TRIDIAGONAL_ELEMENT(center_matrix, row, column, value); \ - \ - row = center_index; \ - column = left_index; \ - value = -coeff1 * arr; /* Left */ \ - UPDATE_TRIDIAGONAL_ELEMENT(center_matrix, row, column, value); \ - \ - /* Remark: Right is not included here due to the symmetry shift */ \ - \ - row = center_index; \ - column = center_index; \ - value = (coeff1 + coeff2) * arr + (coeff3 + coeff4) * att; /* Center: (Left, Right, Bottom, Top) */ \ - UPDATE_TRIDIAGONAL_ELEMENT(center_matrix, row, column, value); \ - \ - /* Fill matrix row of (i-1,j) */ \ - row = left_index; \ - column = center_index; \ - value = -coeff1 * arr; /* Right */ \ - UPDATE_TRIDIAGONAL_ELEMENT(center_matrix, row, column, value); \ - \ - row = left_index; \ - column = left_index; \ - value = coeff1 * arr; /* Center: (Right) */ \ - UPDATE_TRIDIAGONAL_ELEMENT(center_matrix, row, column, value); \ - \ - /* Fill matrix row of (i+1,j) */ \ - /* Nothing to be done */ \ - \ - /* Fill matrix row of (i,j-1) */ \ - row = bottom_index; \ - column = bottom_index; \ - value = coeff3 * att; /* Center: (Top) */ \ - UPDATE_DIAGONAL_ELEMENT(bottom_matrix, row, column, value); \ - \ - /* Fill matrix row of (i,j+1) */ \ - row = top_index; \ - column = top_index; \ - value = coeff4 * att; /* Center: (Bottom) */ \ - UPDATE_DIAGONAL_ELEMENT(top_matrix, row, column, value); \ - } \ - else { \ - /* i_theta % 2 == 0 */ \ - /* ---------------|| */ \ - /* o o o o || */ \ - /* ---------------|| */ \ - /* o x O x || */ \ - /* ---------------|| */ \ - /* o o o o || */ \ - /* ---------------|| */ \ - auto& center_matrix = radial_diagonal_solver[i_theta / 2]; \ - auto& bottom_matrix = radial_tridiagonal_solver[i_theta_M1 / 2]; \ - auto& top_matrix = radial_tridiagonal_solver[i_theta_P1 / 2]; \ - /* Fill matrix row of (i,j) */ \ - row = center_index; \ - column = center_index; \ - value = 0.25 * (h1 + h2) * (k1 + k2) * coeff_beta * fabs(detDF); /* Center: beta_{i,j} */ \ - UPDATE_DIAGONAL_ELEMENT(center_matrix, row, column, value); \ - \ - row = center_index; \ - column = center_index; \ - value = (coeff1 + coeff2) * arr + (coeff3 + coeff4) * att; /* Center: (Left, Right, Bottom, Top) */ \ - UPDATE_DIAGONAL_ELEMENT(center_matrix, row, column, value); \ - \ - /* Fill matrix row of (i,j-1) */ \ - row = bottom_index; \ - column = bottom_index; \ - value = coeff3 * att; /* Center: (Top) */ \ - UPDATE_TRIDIAGONAL_ELEMENT(bottom_matrix, row, column, value); \ - \ - /* Fill matrix row of (i,j+1) */ \ - row = top_index; \ - column = top_index; \ - value = coeff4 * att; /* Center: (Bottom) */ \ - UPDATE_TRIDIAGONAL_ELEMENT(top_matrix, row, column, value); \ - } \ - } \ - /* ------------------------------------------ */ \ - /* Radial Section: Node on the outer boundary */ \ - /* ------------------------------------------ */ \ - else if (i_r == grid.nr() - 1) { \ - assert(!i_r % 2 == 0); \ - \ - double h1 = grid.radialSpacing(i_r - 1); \ - double k1 = grid.angularSpacing(i_theta - 1); \ - double k2 = grid.angularSpacing(i_theta); \ - double coeff1 = 0.5 * (k1 + k2) / h1; \ - \ - int center_index = i_r - numberSmootherCircles; \ - int left_index = i_r - numberSmootherCircles - 1; \ - \ - if (i_theta & 1) { \ - /* i_theta % 2 == 1 */ \ - /* -----------|| */ \ - /* x o x || */ \ - /* -----------|| */ \ - /* o o O || */ \ - /* -----------|| */ \ - /* x o x || */ \ - /* -----------|| */ \ - auto& center_matrix = radial_tridiagonal_solver[i_theta / 2]; \ - \ - /* Fill matrix row of (i,j) */ \ - row = center_index; \ - column = center_index; \ - value = 1.0; \ - UPDATE_TRIDIAGONAL_ELEMENT(center_matrix, row, column, value); \ - \ - /* Fill matrix row of (i-1,j) */ \ - row = left_index; \ - column = left_index; \ - value = coeff1 * arr; /* Center: (Right) */ \ - UPDATE_TRIDIAGONAL_ELEMENT(center_matrix, row, column, value); \ - } \ - else { \ - /* i_theta % 2 == 0 */ \ - /* -----------|| */ \ - /* o o o || */ \ - /* -----------|| */ \ - /* x o X || */ \ - /* -----------|| */ \ - /* o o o || */ \ - /* -----------|| */ \ - auto& center_matrix = radial_diagonal_solver[i_theta / 2]; \ - \ - /* Fill matrix row of (i,j) */ \ - row = center_index; \ - column = center_index; \ - value = 1.0; \ - UPDATE_DIAGONAL_ELEMENT(center_matrix, row, column, value); \ - \ - /* Fill matrix row of (i-1,j) */ \ - row = left_index; \ - column = left_index; \ - value = coeff1 * arr; /* Center: (Right) */ \ - UPDATE_DIAGONAL_ELEMENT(center_matrix, row, column, value); \ - } \ - } \ - } while (0) +void ExtrapolatedSmootherGive::nodeBuildSmootherGive( + int i_r, int i_theta, const PolarGrid& grid, bool DirBC_Interior, MatrixType& inner_boundary_circle_matrix, + std::vector>& circle_diagonal_solver, + std::vector>& radial_diagonal_solver, + std::vector>& circle_tridiagonal_solver, + std::vector>& radial_tridiagonal_solver, double arr, double att, double art, + double detDF, double coeff_beta) +{ + assert(i_r >= 0 && i_r < grid.nr()); + assert(i_theta >= 0 && i_theta < grid.ntheta()); + + const int numberSmootherCircles = grid.numberSmootherCircles(); + const int lengthSmootherRadial = grid.lengthSmootherRadial(); + + assert(numberSmootherCircles >= 3); + assert(lengthSmootherRadial >= 3); + + int ptr, offset; + int row, column, col; + double value, val; + /* ------------------------------------------ */ + /* Node in the interior of the Circle Section */ + /* ------------------------------------------ */ + if (i_r > 0 && i_r < numberSmootherCircles - 1) { + double h1 = grid.radialSpacing(i_r - 1); + double h2 = grid.radialSpacing(i_r); + double k1 = grid.angularSpacing(i_theta - 1); + double k2 = grid.angularSpacing(i_theta); + + double coeff1 = 0.5 * (k1 + k2) / h1; + double coeff2 = 0.5 * (k1 + k2) / h2; + double coeff3 = 0.5 * (h1 + h2) / k1; + double coeff4 = 0.5 * (h1 + h2) / k2; + + int i_theta_M1 = grid.wrapThetaIndex(i_theta - 1); + int i_theta_P1 = grid.wrapThetaIndex(i_theta + 1); + + int center_index = i_theta; + int left_index = i_theta; + int right_index = i_theta; + int bottom_index = i_theta_M1; + int top_index = i_theta_P1; + /* -------------------------- */ + /* Cyclic Tridiagonal Section */ + /* i_r % 2 == 1 */ + if (i_r & 1) { + /* i_theta % 2 == 1 */ + /* | x | o | x | */ + /* | | | | */ + /* | o | O | o | */ + /* | | | | */ + /* | x | o | x | */ + /* or */ + /* i_theta % 2 == 0 */ + /* | o | o | o | */ + /* | | | | */ + /* | x | O | x | */ + /* | | | | */ + /* | o | o | o | */ + + auto& center_matrix = circle_tridiagonal_solver[i_r / 2]; + auto& left_matrix = circle_diagonal_solver[(i_r - 1) / 2]; + auto& right_matrix = circle_diagonal_solver[(i_r + 1) / 2]; + + /* Fill matrix row of (i,j) */ + row = center_index; + column = center_index; + value = 0.25 * (h1 + h2) * (k1 + k2) * coeff_beta * fabs(detDF); /* Center: beta_{i,j} */ + updateTridiagonalElement(center_matrix, row, column, value); + + row = center_index; + column = bottom_index; + value = -coeff3 * att; /* Bottom */ + updateTridiagonalElement(center_matrix, row, column, value); + + row = center_index; + column = top_index; + value = -coeff4 * att; /* Top */ + updateTridiagonalElement(center_matrix, row, column, value); + + row = center_index; + column = center_index; + value = (coeff1 + coeff2) * arr + (coeff3 + coeff4) * att; /* Center: (Left, Right, Bottom, Top) */ + updateTridiagonalElement(center_matrix, row, column, value); + + /* Fill matrix row of (i,j-1) */ + row = bottom_index; + column = center_index; + value = -coeff3 * att; /* Top */ + updateTridiagonalElement(center_matrix, row, column, value); + + row = bottom_index; + column = bottom_index; + value = coeff3 * att; /* Center: (Top) */ + updateTridiagonalElement(center_matrix, row, column, value); + + /* Fill matrix row of (i,j+1) */ + row = top_index; + column = center_index; + value = -coeff4 * att; /* Bottom */ + updateTridiagonalElement(center_matrix, row, column, value); + + row = top_index; + column = top_index; + value = coeff4 * att; /* Center: (Bottom) */ + updateTridiagonalElement(center_matrix, row, column, value); + + if (i_theta & 1) { + /* i_theta % 2 == 1 */ + /* | x | o | x | */ + /* | | | | */ + /* | o | O | o | */ + /* | | | | */ + /* | x | o | x | */ + + /* Fill matrix row of (i-1,j) */ + if (i_r == 1) { + /* Only in the case of AcrossOrigin */ + if (!DirBC_Interior) { + row = left_index; + ptr = getCircleAscIndex(i_r - 1, i_theta); + + const Stencil& LeftStencil = getStencil(i_r - 1, i_theta); + + offset = LeftStencil[StencilPosition::Center]; + col = left_index; + val = +coeff1 * arr; /* Center: (Right) */ + updateCOOCSRMatrixElement(inner_boundary_circle_matrix, ptr, offset, row, col, val); + } + } + else { + row = left_index; + column = left_index; + value = coeff1 * arr; /* Center: (Right) */ + updateDiagonalElement(left_matrix, row, column, value); + } + + /* Fill matrix row of (i+1,j) */ + row = right_index; + column = right_index; + value = coeff2 * arr; /* Center: (Left) */ + updateDiagonalElement(right_matrix, row, column, value); + } + } + /* ---------------- */ + /* Diagonal Section */ + /* i_r % 2 == 0 */ + else { + /* i_theta % 2 == 1 */ + /* | o | x | o | */ + /* | | | | */ + /* | o | O | o | */ + /* | | | | */ + /* | o | x | o | */ + /* or */ + /* i_theta % 2 == 0 */ + /* | o | o | o | */ + /* | | | | */ + /* | o | X | o | */ + /* | | | | */ + /* | o | o | o | */ + + auto& center_matrix = circle_diagonal_solver[i_r / 2]; + auto& left_matrix = circle_tridiagonal_solver[(i_r - 1) / 2]; + auto& right_matrix = circle_tridiagonal_solver[(i_r + 1) / 2]; + + if (i_theta & 1) { /* i_theta % 2 == 1 */ + /* Fill matrix row of (i,j) */ + row = center_index; + column = center_index; + value = 0.25 * (h1 + h2) * (k1 + k2) * coeff_beta * std::fabs(detDF); /* Center: beta_{i,j} */ + updateDiagonalElement(center_matrix, row, column, value); + + row = center_index; + column = center_index; + value = (coeff1 + coeff2) * arr + (coeff3 + coeff4) * att; /* Center: (Left, Right, Bottom, Top) */ + updateDiagonalElement(center_matrix, row, column, value); + } + else { /* i_theta % 2 == 0 */ + /* Fill matrix row of (i,j) */ + row = center_index; + column = center_index; + value = 1.0; + updateDiagonalElement(center_matrix, row, column, value); + + /* Fill matrix row of (i,j-1) */ + row = bottom_index; + column = bottom_index; + value = coeff3 * att; /* Center: (Top) */ + updateDiagonalElement(center_matrix, row, column, value); + + /* Fill matrix row of (i,j+1) */ + row = top_index; + column = top_index; + value = coeff4 * att; /* Center: (Bottom) */ + updateDiagonalElement(center_matrix, row, column, value); + } + /* Fill matrix row of (i-1,j) */ + row = left_index; + column = left_index; + value = coeff1 * arr; /* Center: (Right) */ + updateTridiagonalElement(left_matrix, row, column, value); + + /* Fill matrix row of (i+1,j) */ + row = right_index; + column = right_index; + value = coeff2 * arr; /* Center: (Left) */ + updateTridiagonalElement(right_matrix, row, column, value); + } + } + /* ------------------------------------------ */ + /* Node in the interior of the Radial Section */ + /* ------------------------------------------ */ + else if (i_r > numberSmootherCircles && i_r < grid.nr() - 2) { + double h1 = grid.radialSpacing(i_r - 1); + double h2 = grid.radialSpacing(i_r); + double k1 = grid.angularSpacing(i_theta - 1); + double k2 = grid.angularSpacing(i_theta); + + double coeff1 = 0.5 * (k1 + k2) / h1; + double coeff2 = 0.5 * (k1 + k2) / h2; + double coeff3 = 0.5 * (h1 + h2) / k1; + double coeff4 = 0.5 * (h1 + h2) / k2; + + int i_theta_M1 = grid.wrapThetaIndex(i_theta - 1); + int i_theta_P1 = grid.wrapThetaIndex(i_theta + 1); + + int center_index = i_r - numberSmootherCircles; + int left_index = i_r - numberSmootherCircles - 1; + int right_index = i_r - numberSmootherCircles + 1; + int bottom_index = i_r - numberSmootherCircles; + int top_index = i_r - numberSmootherCircles; + /* ------------------- */ + /* Tridiagonal Section */ + /* i_theta % 2 == 1 */ + if (i_theta & 1) { + /* i_r % 2 == 1 */ + /* ---------- */ + /* x o x */ + /* ---------- */ + /* o O o */ + /* ---------- */ + /* x o x */ + /* ---------- */ + /* or */ + /* i_r % 2 == 0 */ + /* ---------- */ + /* o x o */ + /* ---------- */ + /* o O o */ + /* ---------- */ + /* o x o */ + /* ---------- */ + + auto& center_matrix = radial_tridiagonal_solver[i_theta / 2]; + auto& bottom_matrix = radial_diagonal_solver[i_theta_M1 / 2]; + auto& top_matrix = radial_diagonal_solver[i_theta_P1 / 2]; + + /* Fill matrix row of (i,j) */ + row = center_index; + column = center_index; + value = 0.25 * (h1 + h2) * (k1 + k2) * coeff_beta * fabs(detDF); /* Center: beta_{i,j} */ + updateTridiagonalElement(center_matrix, row, column, value); + + row = center_index; + column = left_index; + value = -coeff1 * arr; /* Left */ + updateTridiagonalElement(center_matrix, row, column, value); + + row = center_index; + column = right_index; + value = -coeff2 * arr; /* Right */ + updateTridiagonalElement(center_matrix, row, column, value); + + row = center_index; + column = center_index; + value = (coeff1 + coeff2) * arr + (coeff3 + coeff4) * att; /* Center: (Left, Right, Bottom, Top) */ + updateTridiagonalElement(center_matrix, row, column, value); + + /* Fill matrix row of (i-1,j) */ + row = left_index; + column = center_index; + value = -coeff1 * arr; /* Right */ + updateTridiagonalElement(center_matrix, row, column, value); + + row = left_index; + column = left_index; + value = coeff1 * arr; /* Center: (Right) */ + updateTridiagonalElement(center_matrix, row, column, value); + + /* Fill matrix row of (i+1,j) */ + row = right_index; + column = center_index; + value = -coeff2 * arr; /* Left */ + updateTridiagonalElement(center_matrix, row, column, value); + + row = right_index; + column = right_index; + value = coeff2 * arr; /* Center: (Left) */ + updateTridiagonalElement(center_matrix, row, column, value); + + if (i_r & 1) { /* i_r % 2 == 1 */ + /* ---------- */ + /* x o x */ + /* ---------- */ + /* o O o */ + /* ---------- */ + /* x o x */ + /* ---------- */ + /* Fill matrix row of (i,j-1) */ + row = bottom_index; + column = bottom_index; + value = coeff3 * att; /* Center: (Top) */ + updateDiagonalElement(bottom_matrix, row, column, value); + + /* Fill matrix row of (i,j+1) */ + row = top_index; + column = top_index; + value = coeff4 * att; /* Center: (Bottom) */ + updateDiagonalElement(top_matrix, row, column, value); + } + } + /* ---------------- */ + /* Diagonal Section */ + /* i_theta % 2 == 0 */ + else { + /* i_r % 2 == 1 */ + /* ---------- */ + /* o o o */ + /* ---------- */ + /* x O x */ + /* ---------- */ + /* o o o */ + /* ---------- */ + /* or */ + /* i_r % 2 == 0 */ + /* ---------- */ + /* o o o */ + /* ---------- */ + /* o X o */ + /* ---------- */ + /* o o o */ + /* ---------- */ + + auto& center_matrix = radial_diagonal_solver[i_theta / 2]; + auto& bottom_matrix = radial_tridiagonal_solver[i_theta_M1 / 2]; + auto& top_matrix = radial_tridiagonal_solver[i_theta_P1 / 2]; + if (i_r & 1) { /* i_r % 2 == 1 */ + /* Fill matrix row of (i,j) */ + row = center_index; + column = center_index; + value = 0.25 * (h1 + h2) * (k1 + k2) * coeff_beta * std::fabs(detDF); /* Center: beta_{i,j} */ + updateDiagonalElement(center_matrix, row, column, value); + + row = center_index; + column = center_index; + value = (coeff1 + coeff2) * arr + (coeff3 + coeff4) * att; /* Center: (Left, Right, Bottom, Top) */ + updateDiagonalElement(center_matrix, row, column, value); + } + else { /* i_r % 2 == 0 */ + /* Fill matrix row of (i,j) */ + row = center_index; + column = center_index; + value = 1.0; + updateDiagonalElement(center_matrix, row, column, value); + + /* Fill matrix row of (i-1,j) */ + row = left_index; + column = left_index; + value = coeff1 * arr; /* Center: (Right) */ + updateDiagonalElement(center_matrix, row, column, value); + + /* Fill matrix row of (i+1,j) */ + row = right_index; + column = right_index; + value = coeff2 * arr; /* Center: (Left) */ + updateDiagonalElement(center_matrix, row, column, value); + } + /* Fill matrix row of (i,j-1) */ + row = bottom_index; + column = bottom_index; + value = coeff3 * att; /* Center: (Top) */ + updateTridiagonalElement(bottom_matrix, row, column, value); + + /* Fill matrix row of (i,j+1) */ + row = top_index; + column = top_index; + value = coeff4 * att; /* Center: (Bottom) */ + updateTridiagonalElement(top_matrix, row, column, value); + } + } + /* ------------------------------------------ */ + /* Circle Section: Node in the inner boundary */ + /* ------------------------------------------ */ + else if (i_r == 0) { + /* ------------------------------------------------ */ + /* Case 1: Dirichlet boundary on the inner boundary */ + /* ------------------------------------------------ */ + if (DirBC_Interior) { + /* Fill result(i,j) */ + double h2 = grid.radialSpacing(i_r); + double k1 = grid.angularSpacing(i_theta - 1); + double k2 = grid.angularSpacing(i_theta); + + double coeff2 = 0.5 * (k1 + k2) / h2; + + int i_theta_M1 = grid.wrapThetaIndex(i_theta - 1); + int i_theta_P1 = grid.wrapThetaIndex(i_theta + 1); + + auto& center_matrix = inner_boundary_circle_matrix; + auto& right_matrix = circle_tridiagonal_solver[(i_r + 1) / 2]; + + int center_index = i_theta; + int right_index = i_theta; + int bottom_index = i_theta_M1; + int top_index = i_theta_P1; + + /* Fill matrix row of (i,j) */ + row = center_index; + ptr = getCircleAscIndex(i_r, i_theta); + + const Stencil& CenterStencil = getStencil(i_r, i_theta); + + offset = CenterStencil[StencilPosition::Center]; + col = center_index; + val = 1.0; + updateCOOCSRMatrixElement(center_matrix, ptr, offset, row, col, val); + + /* Fill matrix row of (i+1,j) */ + row = right_index; + column = right_index; + value = coeff2 * arr; /* Center: (Left) */ + updateTridiagonalElement(right_matrix, row, column, value); + } + else { + /* ------------------------------------------------------------- */ + /* Case 2: Across origin discretization on the interior boundary */ + /* ------------------------------------------------------------- */ + // h1 gets replaced with 2 * R0. + // (i_r-1,i_theta) gets replaced with (i_r, i_theta + (grid.ntheta()/2)). + // Some more adjustments from the changing the 9-point stencil to the artifical 7-point stencil. + double h1 = 2.0 * grid.radius(0); + double h2 = grid.radialSpacing(i_r); + double k1 = grid.angularSpacing(i_theta - 1); + double k2 = grid.angularSpacing(i_theta); + + double coeff1 = 0.5 * (k1 + k2) / h1; + double coeff2 = 0.5 * (k1 + k2) / h2; + double coeff3 = 0.5 * (h1 + h2) / k1; + double coeff4 = 0.5 * (h1 + h2) / k2; + + const int i_theta_M1 = grid.wrapThetaIndex(i_theta - 1); + const int i_theta_P1 = grid.wrapThetaIndex(i_theta + 1); + const int i_theta_AcrossOrigin = grid.wrapThetaIndex(i_theta + (grid.ntheta() / 2)); + + const int center_index = i_theta; + const int left_index = i_theta_AcrossOrigin; + const int right_index = i_theta; + const int bottom_index = i_theta_M1; + const int top_index = i_theta_P1; + + const int center_nz_index = getCircleAscIndex(i_r, i_theta); + const int bottom_nz_index = getCircleAscIndex(i_r, i_theta_M1); + const int top_nz_index = getCircleAscIndex(i_r, i_theta_P1); + const int left_nz_index = getCircleAscIndex(i_r, i_theta_AcrossOrigin); + + int nz_index; + const Stencil& CenterStencil = getStencil(i_r, i_theta); + + if (i_theta & 1) { + /* i_theta % 2 == 1 */ + /* -| x | o | x | */ + /* -| | | | */ + /* -| O | o | o | */ + /* -| | | | */ + /* -| x | o | x | */ + + auto& center_matrix = inner_boundary_circle_matrix; + auto& right_matrix = circle_tridiagonal_solver[(i_r + 1) / 2]; + auto& left_matrix = inner_boundary_circle_matrix; + /* Fill matrix row of (i,j) */ + row = center_index; + ptr = center_nz_index; + + offset = CenterStencil[StencilPosition::Center]; + col = center_index; + val = 0.25 * (h1 + h2) * (k1 + k2) * coeff_beta * std::fabs(detDF); /* beta_{i,j} */ + updateCOOCSRMatrixElement(center_matrix, ptr, offset, row, col, val); + + offset = CenterStencil[StencilPosition::Left]; + col = left_index; + val = -coeff1 * arr; /* Left */ + updateCOOCSRMatrixElement(center_matrix, ptr, offset, row, col, val); + + offset = CenterStencil[StencilPosition::Center]; + col = center_index; + val = (coeff1 + coeff2) * arr + (coeff3 + coeff4) * att; /* Center: (Left, Right, Bottom, Top) */ + updateCOOCSRMatrixElement(center_matrix, ptr, offset, row, col, val); + + /* Fill matrix row of (i-1,j) */ + /* From view the view of the across origin node, */ + /* the directions are roatated by 180 degrees in the stencil! */ + row = left_index; + ptr = left_nz_index; + + const Stencil& LeftStencil = CenterStencil; + + offset = LeftStencil[StencilPosition::Left]; + col = center_index; + val = -coeff1 * arr; /* Right -> Left*/ + updateCOOCSRMatrixElement(left_matrix, ptr, offset, row, col, val); + + offset = LeftStencil[StencilPosition::Center]; + col = left_index; + val = +coeff1 * arr; /* Center: (Right) -> Center: (Left) */ + updateCOOCSRMatrixElement(left_matrix, ptr, offset, row, col, val); + + /* Fill matrix row of (i+1,j) */ + row = right_index; + column = right_index; + value = coeff2 * arr; /* Center: (Left) */ + updateTridiagonalElement(right_matrix, row, column, value); + } + else { + /* i_theta % 2 == 0 */ + /* -| o | o | o | */ + /* -| | | | */ + /* -| X | o | x | */ + /* -| | | | */ + /* -| o | o | o | */ + + auto& center_matrix = inner_boundary_circle_matrix; + auto& right_matrix = circle_tridiagonal_solver[(i_r + 1) / 2]; + auto& left_matrix = inner_boundary_circle_matrix; + /* Fill matrix row of (i,j) */ + row = center_index; + ptr = center_nz_index; + + offset = CenterStencil[StencilPosition::Center]; + col = center_index; + val = 1.0; + updateCOOCSRMatrixElement(center_matrix, ptr, offset, row, col, val); + + /* Fill matrix row of (i,j-1) */ + row = bottom_index; + ptr = bottom_nz_index; + + const Stencil& BottomStencil = CenterStencil; + + offset = BottomStencil[StencilPosition::Center]; + col = bottom_index; + val = +coeff3 * att; /* Center: (Top) */ + updateCOOCSRMatrixElement(center_matrix, ptr, offset, row, col, val); + + /* Fill matrix row of (i,j+1) */ + row = top_index; + ptr = top_nz_index; + + const Stencil& TopStencil = CenterStencil; + + offset = TopStencil[StencilPosition::Center]; + col = top_index; + val = +coeff4 * att; /* Center: (Bottom) */ + updateCOOCSRMatrixElement(center_matrix, ptr, offset, row, col, val); + + /* Fill matrix row of (i+1,j) */ + row = right_index; + column = right_index; + value = coeff2 * arr; /* Center: (Left) */ + updateTridiagonalElement(right_matrix, row, column, value); + } + } + } + /* ------------------------------------------- */ + /* Circle Section: Node next to radial section */ + /* ------------------------------------------- */ + else if (i_r == numberSmootherCircles - 1) { + assert(i_r > 1); + + double h1 = grid.radialSpacing(i_r - 1); + double h2 = grid.radialSpacing(i_r); + double k1 = grid.angularSpacing(i_theta - 1); + double k2 = grid.angularSpacing(i_theta); + + double coeff1 = 0.5 * (k1 + k2) / h1; + double coeff2 = 0.5 * (k1 + k2) / h2; + double coeff3 = 0.5 * (h1 + h2) / k1; + double coeff4 = 0.5 * (h1 + h2) / k2; + + int i_theta_M1 = grid.wrapThetaIndex(i_theta - 1); + int i_theta_P1 = grid.wrapThetaIndex(i_theta + 1); + + int center_index = i_theta; + int left_index = i_theta; + int right_index = 0; + int bottom_index = i_theta_M1; + int top_index = i_theta_P1; + + if (i_r & 1) { + if (i_theta & 1) { + /* i_r % 2 == 1 and i_theta % 2 == 1 */ + /* | o | x | o || x o x o */ + /* | | | || -------------- */ + /* | o | o | O || o o o o */ + /* | | | || -------------- */ + /* | o | x | o || x o x o */ + + auto& center_matrix = circle_tridiagonal_solver[i_r / 2]; + auto& left_matrix = circle_diagonal_solver[(i_r - 1) / 2]; + auto& right_matrix = radial_tridiagonal_solver[i_theta / 2]; + + /* Fill matrix row of (i,j) */ + row = center_index; + column = center_index; + value = 0.25 * (h1 + h2) * (k1 + k2) * coeff_beta * std::fabs(detDF); /* Center: beta_{i,j} */ + updateTridiagonalElement(center_matrix, row, column, value); + + row = center_index; + column = bottom_index; + value = -coeff3 * att; /* Bottom */ + updateTridiagonalElement(center_matrix, row, column, value); + + row = center_index; + column = top_index; + value = -coeff4 * att; /* Top */ + updateTridiagonalElement(center_matrix, row, column, value); + + row = center_index; + column = center_index; + value = (coeff1 + coeff2) * arr + (coeff3 + coeff4) * att; /* Center: (Left, Right, Bottom, Top) */ + updateTridiagonalElement(center_matrix, row, column, value); + + /* Fill matrix row of (i,j-1) */ + row = bottom_index; + column = center_index; + value = -coeff3 * att; /* Top */ + updateTridiagonalElement(center_matrix, row, column, value); + + row = bottom_index; + column = bottom_index; + value = coeff3 * att; /* Center: (Top) */ + updateTridiagonalElement(center_matrix, row, column, value); + + /* Fill matrix row of (i,j+1) */ + row = top_index; + column = center_index; + value = -coeff4 * att; /* Bottom */ + updateTridiagonalElement(center_matrix, row, column, value); + + row = top_index; + column = top_index; + value = coeff4 * att; /* Center: (Bottom) */ + updateTridiagonalElement(center_matrix, row, column, value); + + /* Fill matrix row of (i-1,j) */ + row = left_index; + column = left_index; + value = coeff1 * arr; /* Center: (Right) */ + updateDiagonalElement(left_matrix, row, column, value); + + /* Fill matrix row of (i+1,j) */ + row = right_index; + column = right_index; + value = coeff2 * arr; /* Center: (Left) */ + updateTridiagonalElement(right_matrix, row, column, value); + } + else { + /* i_r % 2 == 1 and i_theta % 2 == 0 */ + /* | o | o | o || o o o o */ + /* | | | || -------------- */ + /* | o | x | O || x o x o */ + /* | | | || -------------- */ + /* | o | o | o || o o o o */ + + auto& center_matrix = circle_tridiagonal_solver[i_r / 2]; + auto& left_matrix = circle_diagonal_solver[(i_r - 1) / 2]; + auto& right_matrix = radial_diagonal_solver[i_theta / 2]; + + /* Fill matrix row of (i,j) */ + row = center_index; + column = center_index; + value = 0.25 * (h1 + h2) * (k1 + k2) * coeff_beta * std::fabs(detDF); /* Center: beta_{i,j} */ + updateTridiagonalElement(center_matrix, row, column, value); + + row = center_index; + column = bottom_index; + value = -coeff3 * att; /* Bottom */ + updateTridiagonalElement(center_matrix, row, column, value); + + row = center_index; + column = top_index; + value = -coeff4 * att; /* Top */ + updateTridiagonalElement(center_matrix, row, column, value); + + row = center_index; + column = center_index; + value = (coeff1 + coeff2) * arr + (coeff3 + coeff4) * att; /* Center: (Left, Right, Bottom, Top) */ + updateTridiagonalElement(center_matrix, row, column, value); + + /* Fill matrix row of (i,j-1) */ + row = bottom_index; + column = center_index; + value = -coeff3 * att; /* Top */ + updateTridiagonalElement(center_matrix, row, column, value); + + row = bottom_index; + column = bottom_index; + value = coeff3 * att; /* Center: (Top) */ + updateTridiagonalElement(center_matrix, row, column, value); + + /* Fill matrix row of (i,j+1) */ + row = top_index; + column = center_index; + value = -coeff4 * att; /* Bottom */ + updateTridiagonalElement(center_matrix, row, column, value); + + row = top_index; + column = top_index; + value = coeff4 * att; /* Center: (Bottom) */ + updateTridiagonalElement(center_matrix, row, column, value); + } + } + else { + if (i_theta & 1) { + /* i_r % 2 == 0 and i_theta % 2 == 1 */ + /* | x | o | x || o x o x */ + /* | | | || -------------- */ + /* | o | o | O || o o o o */ + /* | | | || -------------- */ + /* | x | o | x || o x o x */ + + auto& center_matrix = circle_diagonal_solver[i_r / 2]; + auto& left_matrix = circle_tridiagonal_solver[(i_r - 1) / 2]; + auto& right_matrix = radial_tridiagonal_solver[i_theta / 2]; + + /* Fill matrix row of (i,j) */ + row = center_index; + column = center_index; + value = 0.25 * (h1 + h2) * (k1 + k2) * coeff_beta * std::fabs(detDF); /* Center: beta_{i,j} */ + updateDiagonalElement(center_matrix, row, column, value); + + row = center_index; + column = center_index; + value = (coeff1 + coeff2) * arr + (coeff3 + coeff4) * att; /* Center: (Left, Right, Bottom, Top) */ + updateDiagonalElement(center_matrix, row, column, value); + + /* Fill matrix row of (i-1,j) */ + row = left_index; + column = left_index; + value = coeff1 * arr; /* Center: (Right) */ + updateTridiagonalElement(left_matrix, row, column, value); + + /* Fill matrix row of (i+1,j) */ + row = right_index; + column = right_index; + value = coeff2 * arr; /* Center: (Left) */ + updateTridiagonalElement(right_matrix, row, column, value); + } + else { + /* i_r % 2 == 0 and i_theta % 2 == 0 */ + /* | o | o | o || o o o o */ + /* | | | || -------------- */ + /* | x | o | X || o x o x */ + /* | | | || -------------- */ + /* | o | o | o || o o o o */ + + auto& center_matrix = circle_diagonal_solver[i_r / 2]; + auto& left_matrix = circle_tridiagonal_solver[(i_r - 1) / 2]; + auto& right_matrix = radial_diagonal_solver[i_theta / 2]; + + /* Fill matrix row of (i,j) */ + row = center_index; + column = center_index; + value = 1.0; + updateDiagonalElement(center_matrix, row, column, value); + + /* Fill matrix row of (i,j-1) */ + row = bottom_index; + column = bottom_index; + value = coeff3 * att; /* Center: (Top) */ + updateDiagonalElement(center_matrix, row, column, value); + + /* Fill matrix row of (i,j+1) */ + row = top_index; + column = top_index; + value = coeff4 * att; /* Center: (Bottom) */ + updateDiagonalElement(center_matrix, row, column, value); + + /* Fill matrix row of (i-1,j) */ + row = left_index; + column = left_index; + value = coeff1 * arr; /* Center: (Right) */ + updateTridiagonalElement(left_matrix, row, column, value); + + /* Fill matrix row of (i+1,j) */ + row = right_index; + column = right_index; + value = coeff2 * arr; /* Center: (Left) */ + updateDiagonalElement(right_matrix, row, column, value); + } + } + } + /* --------------------------------------------- */ + /* Radial Section: Node next to circular section */ + /* --------------------------------------------- */ + else if (i_r == numberSmootherCircles) { + double h1 = grid.radialSpacing(i_r - 1); + double h2 = grid.radialSpacing(i_r); + double k1 = grid.angularSpacing(i_theta - 1); + double k2 = grid.angularSpacing(i_theta); + + double coeff1 = 0.5 * (k1 + k2) / h1; + double coeff2 = 0.5 * (k1 + k2) / h2; + double coeff3 = 0.5 * (h1 + h2) / k1; + double coeff4 = 0.5 * (h1 + h2) / k2; + + const int i_theta_M1 = grid.wrapThetaIndex(i_theta - 1); + const int i_theta_P1 = grid.wrapThetaIndex(i_theta + 1); + + const int center_index = i_r - numberSmootherCircles; + const int left_index = i_theta; + const int right_index = i_r - numberSmootherCircles + 1; + const int bottom_index = i_r - numberSmootherCircles; + const int top_index = i_r - numberSmootherCircles; + + if (i_theta & 1) { + if (i_r & 1) { + /* i_theta % 2 == 1 and i_r % 2 == 1 */ + /* | x | o | x || o x o x */ + /* | | | || -------------- */ + /* | o | o | o || O o o o */ + /* | | | || -------------- */ + /* | x | o | x || o x o x */ + + auto& center_matrix = radial_tridiagonal_solver[i_theta / 2]; + auto& bottom_matrix = radial_diagonal_solver[i_theta_M1 / 2]; + auto& top_matrix = radial_diagonal_solver[i_theta_P1 / 2]; + auto& left_matrix = circle_diagonal_solver[(i_r - 1) / 2]; + + /* Fill matrix row of (i,j) */ + row = center_index; + column = center_index; + value = 0.25 * (h1 + h2) * (k1 + k2) * coeff_beta * std::fabs(detDF); /* Center: beta_{i,j} */ + updateTridiagonalElement(center_matrix, row, column, value); + + row = center_index; + column = right_index; + value = -coeff2 * arr; /* Right */ + updateTridiagonalElement(center_matrix, row, column, value); + + row = center_index; + column = center_index; + value = (coeff1 + coeff2) * arr + (coeff3 + coeff4) * att; /* Center: (Left, Right, Bottom, Top) */ + updateTridiagonalElement(center_matrix, row, column, value); + + /* Fill matrix row of (i-1,j) */ + row = left_index; + column = left_index; + value = coeff1 * arr; /* Center: (Right) */ + updateDiagonalElement(left_matrix, row, column, value); + + /* Fill matrix row of (i+1,j) */ + row = right_index; + column = center_index; + value = -coeff2 * arr; /* Left */ + updateTridiagonalElement(center_matrix, row, column, value); + + row = right_index; + column = right_index; + value = coeff2 * arr; /* Center: (Left) */ + updateTridiagonalElement(center_matrix, row, column, value); + + /* Fill matrix row of (i,j-1) */ + row = bottom_index; + column = bottom_index; + value = coeff3 * att; /* Center: (Top) */ + updateDiagonalElement(bottom_matrix, row, column, value); + + /* Fill matrix row of (i,j+1) */ + row = top_index; + column = top_index; + value = coeff4 * att; /* Center: (Bottom) */ + updateDiagonalElement(top_matrix, row, column, value); + } + else { + /* i_theta % 2 == 1 and i_r % 2 == 0 */ + /* | o | x | o || x o x o */ + /* | | | || -------------- */ + /* | o | o | o || O o o o */ + /* | | | || -------------- */ + /* | o | x | o || x o x o */ + + auto& center_matrix = radial_tridiagonal_solver[i_theta / 2]; + auto& bottom_matrix = radial_diagonal_solver[i_theta_M1 / 2]; + auto& top_matrix = radial_diagonal_solver[i_theta_P1 / 2]; + auto& left_matrix = circle_tridiagonal_solver[(i_r - 1) / 2]; + + /* Fill matrix row of (i,j) */ + row = center_index; + column = center_index; + value = 0.25 * (h1 + h2) * (k1 + k2) * coeff_beta * std::fabs(detDF); /* Center: beta_{i,j} */ + updateTridiagonalElement(center_matrix, row, column, value); + + row = center_index; + column = right_index; + value = -coeff2 * arr; /* Right */ + updateTridiagonalElement(center_matrix, row, column, value); + + row = center_index; + column = center_index; + value = (coeff1 + coeff2) * arr + (coeff3 + coeff4) * att; /* Center: (Left, Right, Bottom, Top) */ + updateTridiagonalElement(center_matrix, row, column, value); + + /* Fill matrix row of (i-1,j) */ + row = left_index; + column = left_index; + value = coeff1 * arr; /* Center: (Right) */ + updateTridiagonalElement(left_matrix, row, column, value); + + /* Fill matrix row of (i+1,j) */ + row = right_index; + column = center_index; + value = -coeff2 * arr; /* Left */ + updateTridiagonalElement(center_matrix, row, column, value); + + row = right_index; + column = right_index; + value = coeff2 * arr; /* Center: (Left) */ + updateTridiagonalElement(center_matrix, row, column, value); + } + } + else { + if (i_r & 1) { + /* i_theta % 2 == 0 and i_r % 2 == 1 */ + /* | o | o | o || o o o o */ + /* | | | || -------------- */ + /* | x | o | x || O x o x */ + /* | | | || -------------- */ + /* | o | o | o || o o o o */ + + auto& center_matrix = radial_diagonal_solver[i_theta / 2]; + auto& bottom_matrix = radial_tridiagonal_solver[i_theta_M1 / 2]; + auto& top_matrix = radial_tridiagonal_solver[i_theta_P1 / 2]; + auto& left_matrix = circle_diagonal_solver[(i_r - 1) / 2]; + + /* Fill matrix row of (i,j) */ + row = center_index; + column = center_index; + value = 0.25 * (h1 + h2) * (k1 + k2) * coeff_beta * std::fabs(detDF); /* Center: beta_{i,j} */ + updateDiagonalElement(center_matrix, row, column, value); + + row = center_index; + column = center_index; + value = (coeff1 + coeff2) * arr + (coeff3 + coeff4) * att; /* Center: (Left, Right, Bottom, Top) */ + updateDiagonalElement(center_matrix, row, column, value); + + /* Fill matrix row of (i,j-1) */ + row = bottom_index; + column = bottom_index; + value = coeff3 * att; /* Center: (Top) */ + updateTridiagonalElement(bottom_matrix, row, column, value); + + /* Fill matrix row of (i,j+1) */ + row = top_index; + column = top_index; + value = coeff4 * att; /* Center: (Bottom) */ + updateTridiagonalElement(top_matrix, row, column, value); + } + else { + /* i_theta % 2 == 0 and i_r % 2 == 0 */ + /* | o | o | o || o o o o */ + /* | | | || -------------- */ + /* | o | x | o || X o x o */ + /* | | | || -------------- */ + /* | o | o | o || o o o o */ + + auto& center_matrix = radial_diagonal_solver[i_theta / 2]; + auto& bottom_matrix = radial_tridiagonal_solver[i_theta_M1 / 2]; + auto& top_matrix = radial_tridiagonal_solver[i_theta_P1 / 2]; + auto& left_matrix = circle_tridiagonal_solver[(i_r - 1) / 2]; + + /* Fill matrix row of (i,j) */ + + row = center_index; + column = center_index; + value = 1.0; + updateDiagonalElement(center_matrix, row, column, value); + + /* Fill matrix row of (i-1,j) */ + row = left_index; + column = left_index; + value = coeff1 * arr; /* Center: (Right) */ + updateTridiagonalElement(left_matrix, row, column, value); + + /* Fill matrix row of (i+1,j) */ + row = right_index; + column = right_index; + value = coeff2 * arr; /* Center: (Left) */ + updateDiagonalElement(center_matrix, row, column, value); + + /* Fill matrix row of (i,j-1) */ + row = bottom_index; + column = bottom_index; + value = coeff3 * att; /* Center: (Top) */ + updateTridiagonalElement(bottom_matrix, row, column, value); + + /* Fill matrix row of (i,j+1) */ + row = top_index; + column = top_index; + value = coeff4 * att; /* Center: (Bottom) */ + updateTridiagonalElement(top_matrix, row, column, value); + } + } + } + /* ------------------------------------------- */ + /* Radial Section: Node next to outer boundary */ + /* ------------------------------------------- */ + else if (i_r == grid.nr() - 2) { + assert(i_r % 2 == 1); + + double h1 = grid.radialSpacing(i_r - 1); + double h2 = grid.radialSpacing(i_r); + double k1 = grid.angularSpacing(i_theta - 1); + double k2 = grid.angularSpacing(i_theta); + + double coeff1 = 0.5 * (k1 + k2) / h1; + double coeff2 = 0.5 * (k1 + k2) / h2; + double coeff3 = 0.5 * (h1 + h2) / k1; + double coeff4 = 0.5 * (h1 + h2) / k2; + + int i_theta_M1 = grid.wrapThetaIndex(i_theta - 1); + int i_theta_P1 = grid.wrapThetaIndex(i_theta + 1); + + int center_index = i_r - numberSmootherCircles; + int left_index = i_r - numberSmootherCircles - 1; + int right_index = i_r - numberSmootherCircles + 1; + int bottom_index = i_r - numberSmootherCircles; + int top_index = i_r - numberSmootherCircles; + + if (i_theta & 1) { + /* i_theta % 2 == 1 */ + /* ---------------|| */ + /* o x o x || */ + /* ---------------|| */ + /* o o O o || */ + /* ---------------|| */ + /* o x o x || */ + /* ---------------|| */ + auto& center_matrix = radial_tridiagonal_solver[i_theta / 2]; + auto& bottom_matrix = radial_diagonal_solver[i_theta_M1 / 2]; + auto& top_matrix = radial_diagonal_solver[i_theta_P1 / 2]; + /* Fill matrix row of (i,j) */ + row = center_index; + column = center_index; + value = 0.25 * (h1 + h2) * (k1 + k2) * coeff_beta * std::fabs(detDF); /* Center: beta_{i,j} */ + updateTridiagonalElement(center_matrix, row, column, value); + + row = center_index; + column = left_index; + value = -coeff1 * arr; /* Left */ + updateTridiagonalElement(center_matrix, row, column, value); + + /* Remark: Right is not included here due to the symmetry shift */ + + row = center_index; + column = center_index; + value = (coeff1 + coeff2) * arr + (coeff3 + coeff4) * att; /* Center: (Left, Right, Bottom, Top) */ + updateTridiagonalElement(center_matrix, row, column, value); + + /* Fill matrix row of (i-1,j) */ + row = left_index; + column = center_index; + value = -coeff1 * arr; /* Right */ + updateTridiagonalElement(center_matrix, row, column, value); + + row = left_index; + column = left_index; + value = coeff1 * arr; /* Center: (Right) */ + updateTridiagonalElement(center_matrix, row, column, value); + + /* Fill matrix row of (i+1,j) */ + /* Nothing to be done */ + + /* Fill matrix row of (i,j-1) */ + row = bottom_index; + column = bottom_index; + value = coeff3 * att; /* Center: (Top) */ + updateDiagonalElement(bottom_matrix, row, column, value); + + /* Fill matrix row of (i,j+1) */ + row = top_index; + column = top_index; + value = coeff4 * att; /* Center: (Bottom) */ + updateDiagonalElement(top_matrix, row, column, value); + } + else { + /* i_theta % 2 == 0 */ + /* ---------------|| */ + /* o o o o || */ + /* ---------------|| */ + /* o x O x || */ + /* ---------------|| */ + /* o o o o || */ + /* ---------------|| */ + auto& center_matrix = radial_diagonal_solver[i_theta / 2]; + auto& bottom_matrix = radial_tridiagonal_solver[i_theta_M1 / 2]; + auto& top_matrix = radial_tridiagonal_solver[i_theta_P1 / 2]; + /* Fill matrix row of (i,j) */ + row = center_index; + column = center_index; + value = 0.25 * (h1 + h2) * (k1 + k2) * coeff_beta * std::fabs(detDF); /* Center: beta_{i,j} */ + updateDiagonalElement(center_matrix, row, column, value); + + row = center_index; + column = center_index; + value = (coeff1 + coeff2) * arr + (coeff3 + coeff4) * att; /* Center: (Left, Right, Bottom, Top) */ + updateDiagonalElement(center_matrix, row, column, value); + + /* Fill matrix row of (i,j-1) */ + row = bottom_index; + column = bottom_index; + value = coeff3 * att; /* Center: (Top) */ + updateTridiagonalElement(bottom_matrix, row, column, value); + + /* Fill matrix row of (i,j+1) */ + row = top_index; + column = top_index; + value = coeff4 * att; /* Center: (Bottom) */ + updateTridiagonalElement(top_matrix, row, column, value); + } + } + /* ------------------------------------------ */ + /* Radial Section: Node on the outer boundary */ + /* ------------------------------------------ */ + else if (i_r == grid.nr() - 1) { + assert(!i_r % 2 == 0); + + double h1 = grid.radialSpacing(i_r - 1); + double k1 = grid.angularSpacing(i_theta - 1); + double k2 = grid.angularSpacing(i_theta); + + double coeff1 = 0.5 * (k1 + k2) / h1; + + int center_index = i_r - numberSmootherCircles; + int left_index = i_r - numberSmootherCircles - 1; + + if (i_theta & 1) { + /* i_theta % 2 == 1 */ + /* -----------|| */ + /* x o x || */ + /* -----------|| */ + /* o o O || */ + /* -----------|| */ + /* x o x || */ + /* -----------|| */ + auto& center_matrix = radial_tridiagonal_solver[i_theta / 2]; + + /* Fill matrix row of (i,j) */ + row = center_index; + column = center_index; + value = 1.0; + updateTridiagonalElement(center_matrix, row, column, value); + + /* Fill matrix row of (i-1,j) */ + row = left_index; + column = left_index; + value = coeff1 * arr; /* Center: (Right) */ + updateTridiagonalElement(center_matrix, row, column, value); + } + else { + /* i_theta % 2 == 0 */ + /* -----------|| */ + /* o o o || */ + /* -----------|| */ + /* x o X || */ + /* -----------|| */ + /* o o o || */ + /* -----------|| */ + auto& center_matrix = radial_diagonal_solver[i_theta / 2]; + + /* Fill matrix row of (i,j) */ + row = center_index; + column = center_index; + value = 1.0; + updateDiagonalElement(center_matrix, row, column, value); + + /* Fill matrix row of (i-1,j) */ + row = left_index; + column = left_index; + value = coeff1 * arr; /* Center: (Right) */ + updateDiagonalElement(center_matrix, row, column, value); + } + } +} void ExtrapolatedSmootherGive::buildAscCircleSection(const int i_r) { @@ -1237,9 +1243,9 @@ void ExtrapolatedSmootherGive::buildAscCircleSection(const int i_r) level_cache_.obtainValues(i_r, i_theta, global_index, r, theta, coeff_beta, arr, att, art, detDF); // Build Asc at the current node - NODE_BUILD_SMOOTHER_GIVE(i_r, i_theta, grid_, DirBC_Interior_, inner_boundary_circle_matrix_, - circle_diagonal_solver_, circle_tridiagonal_solver_, radial_diagonal_solver_, - radial_tridiagonal_solver_); + nodeBuildSmootherGive(i_r, i_theta, grid_, DirBC_Interior_, inner_boundary_circle_matrix_, + circle_diagonal_solver_, radial_diagonal_solver_, circle_tridiagonal_solver_, + radial_tridiagonal_solver_, arr, att, art, detDF, coeff_beta); } } @@ -1254,9 +1260,9 @@ void ExtrapolatedSmootherGive::buildAscRadialSection(const int i_theta) level_cache_.obtainValues(i_r, i_theta, global_index, r, theta, coeff_beta, arr, att, art, detDF); // Build Asc at the current node - NODE_BUILD_SMOOTHER_GIVE(i_r, i_theta, grid_, DirBC_Interior_, inner_boundary_circle_matrix_, - circle_diagonal_solver_, circle_tridiagonal_solver_, radial_diagonal_solver_, - radial_tridiagonal_solver_); + nodeBuildSmootherGive(i_r, i_theta, grid_, DirBC_Interior_, inner_boundary_circle_matrix_, + circle_diagonal_solver_, radial_diagonal_solver_, circle_tridiagonal_solver_, + radial_tridiagonal_solver_, arr, att, art, detDF, coeff_beta); } } diff --git a/src/ExtrapolatedSmoother/ExtrapolatedSmootherGive/smootherSolver.cpp b/src/ExtrapolatedSmoother/ExtrapolatedSmootherGive/smootherSolver.cpp index 8124c66f..d960640b 100644 --- a/src/ExtrapolatedSmoother/ExtrapolatedSmootherGive/smootherSolver.cpp +++ b/src/ExtrapolatedSmoother/ExtrapolatedSmootherGive/smootherSolver.cpp @@ -2,817 +2,888 @@ #include "../../../include/common/geometry_helper.h" -// The current position is marked with a ~ symbol. -#define NODE_APPLY_ASC_ORTHO_CIRCLE_GIVE(i_r, i_theta, r, theta, grid, DirBC_Interior, smoother_color, x, rhs, temp, \ - arr, att, art, detDF, coeff_beta) \ - do { \ - assert(i_r >= 0 && i_r <= grid_.numberSmootherCircles()); \ - bool isOddNumberSmootherCircles = (grid.numberSmootherCircles() & 1); \ - bool isOddRadialIndex = (i_r & 1); \ - SmootherColor node_color = \ - (isOddNumberSmootherCircles == isOddRadialIndex) ? SmootherColor::White : SmootherColor::Black; \ - /* -------------------- */ \ - /* Node in the interior */ \ - /* -------------------- */ \ - if (i_r > 0 && i_r < grid.numberSmootherCircles()) { \ - double h1 = grid.radialSpacing(i_r - 1); \ - double h2 = grid.radialSpacing(i_r); \ - double k1 = grid.angularSpacing(i_theta - 1); \ - double k2 = grid.angularSpacing(i_theta); \ - double coeff1 = 0.5 * (k1 + k2) / h1; \ - double coeff2 = 0.5 * (k1 + k2) / h2; \ - double coeff3 = 0.5 * (h1 + h2) / k1; \ - double coeff4 = 0.5 * (h1 + h2) / k2; \ - /* -------------------- */ \ - /* Inside Section Parts */ \ - if (node_color == smoother_color) { \ - if (i_r & 1) { \ - if (i_theta & 1) { \ - /* i_r % 2 == 1 and i_theta % 2 == 1 */ \ - /* | x | o | x | */ \ - /* | | | | */ \ - /* | o | O | o | */ \ - /* | | | | */ \ - /* | x | o | x | */ \ - \ - /* Fill temp(i,j) */ \ - temp[grid.index(i_r, i_theta)] -= \ - (-coeff1 * arr * x[grid.index(i_r - 1, i_theta)] /* Left */ \ - - coeff2 * arr * x[grid.index(i_r + 1, i_theta)] /* Right */ \ - ); \ - /* Fill temp(i,j-1) */ \ - temp[grid.index(i_r, i_theta - 1)] -= \ - (-0.25 * art * x[grid.index(i_r + 1, i_theta)] /* Top Right */ \ - + 0.25 * art * x[grid.index(i_r - 1, i_theta)]); /* Top Left */ \ - /* Fill temp(i,j+1) */ \ - temp[grid.index(i_r, i_theta + 1)] -= \ - (+0.25 * art * x[grid.index(i_r + 1, i_theta)] /* Bottom Right */ \ - - 0.25 * art * x[grid.index(i_r - 1, i_theta)]); /* Bottom Left */ \ - } \ - else { \ - /* i_r % 2 == 1 and i_theta % 2 == 0 */ \ - /* | o | o | o | */ \ - /* | | | | */ \ - /* | x | O | x | */ \ - /* | | | | */ \ - /* | o | o | o | */ \ - /* Fill temp(i,j) */ \ - temp[grid.index(i_r, i_theta)] -= \ - (-coeff1 * arr * x[grid.index(i_r - 1, i_theta)] /* Left */ \ - - coeff2 * arr * x[grid.index(i_r + 1, i_theta)] /* Right */ \ - ); \ - /* Fill temp(i,j-1) */ \ - temp[grid.index(i_r, i_theta - 1)] -= \ - (-0.25 * art * x[grid.index(i_r + 1, i_theta)] /* Top Right */ \ - + 0.25 * art * x[grid.index(i_r - 1, i_theta)]); /* Top Left */ \ - /* Fill temp(i,j+1) */ \ - temp[grid.index(i_r, i_theta + 1)] -= \ - (+0.25 * art * x[grid.index(i_r + 1, i_theta)] /* Bottom Right */ \ - - 0.25 * art * x[grid.index(i_r - 1, i_theta)]); /* Bottom Left */ \ - } \ - } \ - else { \ - if (i_theta & 1) { \ - /* i_r % 2 == 0 and i_theta % 2 == 1 */ \ - /* | o | x | o | */ \ - /* | | | | */ \ - /* | o | O | o | */ \ - /* | | | | */ \ - /* | o | x | o | */ \ - /* Fill temp(i,j) */ \ - temp[grid.index(i_r, i_theta)] -= \ - (-coeff1 * arr * x[grid.index(i_r - 1, i_theta)] /* Left */ \ - - coeff2 * arr * x[grid.index(i_r + 1, i_theta)] /* Right */ \ - - coeff3 * att * x[grid.index(i_r, i_theta - 1)] /* Bottom */ \ - - coeff4 * att * x[grid.index(i_r, i_theta + 1)] /* Top */ \ - ); \ - } \ - else { \ - /* i_r % 2 == 0 and i_theta % 2 == 0 */ \ - /* | o | o | o | */ \ - /* | | | | */ \ - /* | o | X | o | */ \ - /* | | | | */ \ - /* | o | o | o | */ \ - \ - /* Fill temp(i,j-1) */ \ - temp[grid.index(i_r, i_theta - 1)] -= \ - (-coeff3 * att * x[grid.index(i_r, i_theta)] /* Top */ \ - - 0.25 * art * x[grid.index(i_r + 1, i_theta)] /* Top Right */ \ - + 0.25 * art * x[grid.index(i_r - 1, i_theta)]); /* Top Left */ \ - /* Fill temp(i,j+1) */ \ - temp[grid.index(i_r, i_theta + 1)] -= \ - (-coeff4 * att * x[grid.index(i_r, i_theta)] /* Bottom */ \ - + 0.25 * art * x[grid.index(i_r + 1, i_theta)] /* Bottom Right */ \ - - 0.25 * art * x[grid.index(i_r - 1, i_theta)]); /* Bottom Left */ \ - } \ - } \ - } \ - /* --------------------- */ \ - /* Outside Section Parts */ \ - else if (node_color != smoother_color) { \ - if (i_r & 1) { \ - if (i_theta & 1) { \ - /* i_r % 2 == 1 and i_theta % 2 == 1 */ \ - /* | x | o | x | */ \ - /* | | | | */ \ - /* | o | O | o | */ \ - /* | | | | */ \ - /* | x | o | x | */ \ - \ - /* Fill temp(i-1,j) */ \ - if (i_r > 1 || !DirBC_Interior) { \ - temp[grid.index(i_r - 1, i_theta)] -= \ - (-coeff1 * arr * x[grid.index(i_r, i_theta)] /* Right */ \ - - 0.25 * art * x[grid.index(i_r, i_theta + 1)] /* Top Right */ \ - + 0.25 * art * x[grid.index(i_r, i_theta - 1)]); /* Bottom Right */ \ - } \ - /* Fill temp(i+1,j) */ \ - if (i_r < grid.numberSmootherCircles() - 1) { \ - temp[grid.index(i_r + 1, i_theta)] -= \ - (-coeff2 * arr * x[grid.index(i_r, i_theta)] /* Left */ \ - + 0.25 * art * x[grid.index(i_r, i_theta + 1)] /* Top Left */ \ - - 0.25 * art * x[grid.index(i_r, i_theta - 1)]); /* Bottom Left */ \ - } \ - } \ - else { \ - /* i_r % 2 == 1 and i_theta % 2 == 0 */ \ - /* | o | o | o | */ \ - /* | | | | */ \ - /* | x | O | x | */ \ - /* | | | | */ \ - /* | o | o | o | */ \ - \ - /* Nothing to do! */ \ - } \ - } \ - else { \ - if (i_theta & 1) { \ - /* i_r % 2 == 0 and i_theta % 2 == 1 */ \ - /* | o | x | o | */ \ - /* | | | | */ \ - /* | o | O | o | */ \ - /* | | | | */ \ - /* | o | x | o | */ \ - \ - /* Fill temp(i-1,j) */ \ - if (i_r > 1 || !DirBC_Interior) { \ - temp[grid.index(i_r - 1, i_theta)] -= \ - (-coeff1 * arr * x[grid.index(i_r, i_theta)] /* Right */ \ - - 0.25 * art * x[grid.index(i_r, i_theta + 1)] /* Top Right */ \ - + 0.25 * art * x[grid.index(i_r, i_theta - 1)]); /* Bottom Right */ \ - } \ - /* Fill temp(i+1,j) */ \ - if (i_r < grid.numberSmootherCircles() - 1) { \ - temp[grid.index(i_r + 1, i_theta)] -= \ - (-coeff2 * arr * x[grid.index(i_r, i_theta)] /* Left */ \ - + 0.25 * art * x[grid.index(i_r, i_theta + 1)] /* Top Left */ \ - - 0.25 * art * x[grid.index(i_r, i_theta - 1)]); /* Bottom Left */ \ - } \ - } \ - else { \ - /* i_r % 2 == 0 and i_theta % 2 == 0 */ \ - /* | o | o | o | */ \ - /* | | | | */ \ - /* | o | X | o | */ \ - /* | | | | */ \ - /* | o | o | o | */ \ - \ - /* Fill temp(i-1,j) */ \ - if (i_r > 1 || !DirBC_Interior) { \ - temp[grid.index(i_r - 1, i_theta)] -= \ - (-coeff1 * arr * x[grid.index(i_r, i_theta)] /* Right */ \ - - 0.25 * art * x[grid.index(i_r, i_theta + 1)] /* Top Right */ \ - + 0.25 * art * x[grid.index(i_r, i_theta - 1)]); /* Bottom Right */ \ - } \ - /* Fill temp(i+1,j) */ \ - if (i_r < grid.numberSmootherCircles() - 1) { \ - temp[grid.index(i_r + 1, i_theta)] -= \ - (-coeff2 * arr * x[grid.index(i_r, i_theta)] /* Left */ \ - + 0.25 * art * x[grid.index(i_r, i_theta + 1)] /* Top Left */ \ - - 0.25 * art * x[grid.index(i_r, i_theta - 1)]); /* Bottom Left */ \ - } \ - } \ - } \ - } \ - } \ - /* -------------------- */ \ - /* Node on the boundary */ \ - /* -------------------- */ \ - else if (i_r == 0) { \ - /* ------------------------------------------------ */ \ - /* Case 1: Dirichlet boundary on the inner boundary */ \ - /* ------------------------------------------------ */ \ - if (DirBC_Interior) { \ - double h2 = grid.radialSpacing(i_r); \ - double k1 = grid.angularSpacing(i_theta - 1); \ - double k2 = grid.angularSpacing(i_theta); \ - double coeff2 = 0.5 * (k1 + k2) / h2; \ - /* -------------------- */ \ - /* Inside Section Parts */ \ - if (node_color == smoother_color) { \ - if (i_theta & 1) { \ - /* i_theta % 2 == 1 */ \ - /* || x | o | x | */ \ - /* || | | | */ \ - /* || O | o | o | */ \ - /* || | | | */ \ - /* || x | o | x | */ \ - /* Nothing to do! */ \ - } \ - else { \ - /* i_theta % 2 == 0 */ \ - /* || o | o | o | */ \ - /* || | | | */ \ - /* || X | o | x | */ \ - /* || | | | */ \ - /* || o | o | o | */ \ - /* Nothing to do! */ \ - } \ - } \ - /* --------------------- */ \ - /* Outside Section Parts */ \ - else if (node_color != smoother_color) { \ - /* Fill temp(i+1,j) */ \ - temp[grid.index(i_r + 1, i_theta)] -= \ - (-coeff2 * arr * x[grid.index(i_r, i_theta)] /* Left */ \ - + 0.25 * art * x[grid.index(i_r, i_theta + 1)] /* Top Left */ \ - - 0.25 * art * x[grid.index(i_r, i_theta - 1)]); /* Bottom Left */ \ - } \ - } \ - else { \ - /* ------------------------------------------------------------- */ \ - /* Case 2: Across origin discretization on the interior boundary */ \ - /* ------------------------------------------------------------- */ \ - /* h1 gets replaced with 2 * R0. */ \ - /* (i_r-1,i_theta) gets replaced with (i_r, i_theta + (grid.ntheta()/2)). */ \ - /* Some more adjustments from the changing the 9-point stencil to the artifical 7-point stencil. */ \ - double h1 = 2.0 * grid.radius(0); \ - double h2 = grid.radialSpacing(i_r); \ - double k1 = grid.angularSpacing(i_theta - 1); \ - double k2 = grid.angularSpacing(i_theta); \ - double coeff1 = 0.5 * (k1 + k2) / h1; \ - double coeff2 = 0.5 * (k1 + k2) / h2; \ - double coeff3 = 0.5 * (h1 + h2) / k1; \ - double coeff4 = 0.5 * (h1 + h2) / k2; \ - /* -------------------- */ \ - /* Inside Section Parts */ \ - if (node_color == smoother_color) { \ - if (i_theta & 1) { \ - /* i_theta % 2 == 1 */ \ - /* -| x | o | x | */ \ - /* -| | | | */ \ - /* -| O | o | o | */ \ - /* -| | | | */ \ - /* -| x | o | x | */ \ - /* Fill temp(i,j) */ \ - temp[grid.index(i_r, i_theta)] -= \ - (/* - coeff1 * arr * x[grid.index(i_r, i_theta + (grid.ntheta()/2))] // Left: Not in Asc_ortho */ \ - -coeff2 * arr * x[grid.index(i_r + 1, i_theta)] /* Right */ \ - - coeff3 * att * x[grid.index(i_r, i_theta - 1)] /* Bottom */ \ - - coeff4 * att * x[grid.index(i_r, i_theta + 1)] /* Top */ \ - ); \ - } \ - else { \ - /* i_theta % 2 == 0 */ \ - /* -| o | o | o | */ \ - /* -| | | | */ \ - /* -| X | o | x | */ \ - /* -| | | | */ \ - /* -| o | o | o | */ \ - \ - /* Fill temp(i,j-1) */ \ - temp[grid.index(i_r, i_theta - 1)] -= \ - (-coeff3 * att * x[grid.index(i_r, i_theta)] /* Top */ \ - - 0.25 * art * x[grid.index(i_r + 1, i_theta)]); /* Top Right */ \ - /* + 0.25 * art * x[grid.index(i_r-1,i_theta)]; // Top Left: REMOVED DUE TO ARTIFICAL 7 POINT STENCIL */ \ - /* Fill temp(i,j+1) */ \ - temp[grid.index(i_r, i_theta + 1)] -= \ - (-coeff4 * att * x[grid.index(i_r, i_theta)] /* Bottom */ \ - + 0.25 * art * x[grid.index(i_r + 1, i_theta)]); /* Bottom Right */ \ - /* - 0.25 * art * x[grid.index(i_r-1,i_theta)]; // Bottom Left: REMOVED DUE TO ARTIFICAL 7 POINT STENCIL */ \ - } \ - } \ - /* --------------------- */ \ - /* Outside Section Parts */ \ - else if (node_color != smoother_color) { \ - /* Fill temp(i+1,j) */ \ - temp[grid.index(i_r + 1, i_theta)] -= \ - (-coeff2 * arr * x[grid.index(i_r, i_theta)] /* Left */ \ - + 0.25 * art * x[grid.index(i_r, i_theta + 1)] /* Top Left */ \ - - 0.25 * art * x[grid.index(i_r, i_theta - 1)]); /* Bottom Left */ \ - } \ - } \ - } \ - /* ----------------------------- */ \ - /* Node next to circular section */ \ - /* ----------------------------- */ \ - else if (i_r == grid.numberSmootherCircles()) { \ - assert(node_color == SmootherColor::White); \ - if (smoother_color == SmootherColor::Black) { \ - double h1 = grid.radialSpacing(i_r - 1); \ - double h2 = grid.radialSpacing(i_r); \ - double k1 = grid.angularSpacing(i_theta - 1); \ - double k2 = grid.angularSpacing(i_theta); \ - double coeff1 = 0.5 * (k1 + k2) / h1; \ - double coeff2 = 0.5 * (k1 + k2) / h2; \ - double coeff3 = 0.5 * (h1 + h2) / k1; \ - double coeff4 = 0.5 * (h1 + h2) / k2; \ - \ - /* i_theta % 2 == 1 and i_r % 2 == 1 */ \ - /* | x | o | x || o x o x */ \ - /* | | | || -------------- */ \ - /* | o | o | o || O o o o */ \ - /* | | | || -------------- */ \ - /* | x | o | x || o x o x */ \ - /* -> Give Left */ \ - \ - /* i_theta % 2 == 1 and i_r % 2 == 0 */ \ - /* | o | x | o || x o x o */ \ - /* | | | || -------------- */ \ - /* | o | o | o || O o o o */ \ - /* | | | || -------------- */ \ - /* | o | x | o || x o x o */ \ - /* -> Give Left */ \ - \ - /* i_theta % 2 == 0 and i_r % 2 == 1 */ \ - /* | o | o | o || o o o o */ \ - /* | | | || -------------- */ \ - /* | x | o | x || O x o x */ \ - /* | | | || -------------- */ \ - /* | o | o | o || o o o o */ \ - /* -> Don't give to the Left! */ \ - \ - /* i_theta % 2 == 0 and i_r % 2 == 0 */ \ - /* | o | o | o || o o o o */ \ - /* | | | || -------------- */ \ - /* | o | x | o || X o x o */ \ - /* | | | || -------------- */ \ - /* | o | o | o || o o o o */ \ - /* -> Give Left */ \ - \ - if (i_theta & 1 || !(i_r & 1)) { \ - /* --------------------- */ \ - /* Outside Section Parts */ \ - /* Fill temp(i-1,j) */ \ - temp[grid.index(i_r - 1, i_theta)] -= \ - (-coeff1 * arr * x[grid.index(i_r, i_theta)] /* Right */ \ - - 0.25 * art * x[grid.index(i_r, i_theta + 1)] /* Top Right */ \ - + 0.25 * art * x[grid.index(i_r, i_theta - 1)]); /* Bottom Right */ \ - } \ - } \ - } \ - } while (0) - -#define NODE_APPLY_ASC_ORTHO_RADIAL_GIVE(i_r, i_theta, r, theta, grid, DirBC_Interior, smoother_color, x, rhs, temp, \ - arr, att, art, detDF, coeff_beta) \ - do { \ - assert(i_r >= grid.numberSmootherCircles() - 1 && i_r < grid.nr()); \ - SmootherColor node_color = (i_theta & 1) ? SmootherColor::White : SmootherColor::Black; \ - /* -------------------- */ \ - /* Node in the interior */ \ - /* -------------------- */ \ - if (i_r > grid.numberSmootherCircles() && i_r < grid.nr() - 2) { \ - double h1 = grid.radialSpacing(i_r - 1); \ - double h2 = grid.radialSpacing(i_r); \ - double k1 = grid.angularSpacing(i_theta - 1); \ - double k2 = grid.angularSpacing(i_theta); \ - double coeff1 = 0.5 * (k1 + k2) / h1; \ - double coeff2 = 0.5 * (k1 + k2) / h2; \ - double coeff3 = 0.5 * (h1 + h2) / k1; \ - double coeff4 = 0.5 * (h1 + h2) / k2; \ - /* -------------------- */ \ - /* Inside Section Parts */ \ - if (node_color == smoother_color) { \ - if (i_theta & 1) { \ - if (i_r & 1) { \ - /* i_theta % 2 == 1 and i_r % 2 == 1 */ \ - /* ---------- */ \ - /* x o x */ \ - /* ---------- */ \ - /* o O o */ \ - /* ---------- */ \ - /* x o x */ \ - /* ---------- */ \ - /* Fill temp(i,j) */ \ - temp[grid.index(i_r, i_theta)] -= \ - (-coeff3 * att * x[grid.index(i_r, i_theta - 1)] /* Bottom */ \ - - coeff4 * att * x[grid.index(i_r, i_theta + 1)] /* Top */ \ - ); \ - /* Fill temp(i-1,j) */ \ - temp[grid.index(i_r - 1, i_theta)] -= \ - (-0.25 * art * x[grid.index(i_r, i_theta + 1)] /* Top Right */ \ - + 0.25 * art * x[grid.index(i_r, i_theta - 1)]); /* Bottom Right */ \ - /* Fill temp(i+1,j) */ \ - temp[grid.index(i_r + 1, i_theta)] -= \ - (+0.25 * art * x[grid.index(i_r, i_theta + 1)] /* Top Left */ \ - - 0.25 * art * x[grid.index(i_r, i_theta - 1)]); /* Bottom Left */ \ - } \ - else { \ - /* i_theta % 2 == 1 and i_r % 2 == 0 */ \ - /* ---------- */ \ - /* o x o */ \ - /* ---------- */ \ - /* o O o */ \ - /* ---------- */ \ - /* o x o */ \ - /* ---------- */ \ - /* Fill temp(i,j) */ \ - temp[grid.index(i_r, i_theta)] -= \ - (-coeff3 * att * x[grid.index(i_r, i_theta - 1)] /* Bottom */ \ - - coeff4 * att * x[grid.index(i_r, i_theta + 1)] /* Top */ \ - ); \ - /* Fill temp(i-1,j) */ \ - temp[grid.index(i_r - 1, i_theta)] -= \ - (-0.25 * art * x[grid.index(i_r, i_theta + 1)] /* Top Right */ \ - + 0.25 * art * x[grid.index(i_r, i_theta - 1)]); /* Bottom Right */ \ - /* Fill temp(i+1,j) */ \ - temp[grid.index(i_r + 1, i_theta)] -= \ - (+0.25 * art * x[grid.index(i_r, i_theta + 1)] /* Top Left */ \ - - 0.25 * art * x[grid.index(i_r, i_theta - 1)]); /* Bottom Left */ \ - } \ - } \ - else { \ - if (i_r & 1) { \ - /* i_theta % 2 == 0 and i_r % 2 == 1 */ \ - /* ---------- */ \ - /* o o o */ \ - /* ---------- */ \ - /* x O x */ \ - /* ---------- */ \ - /* o o o */ \ - /* ---------- */ \ - /* Fill temp(i,j) */ \ - temp[grid.index(i_r, i_theta)] -= \ - (-coeff1 * arr * x[grid.index(i_r - 1, i_theta)] /* Left */ \ - - coeff2 * arr * x[grid.index(i_r + 1, i_theta)] /* Right */ \ - - coeff3 * att * x[grid.index(i_r, i_theta - 1)] /* Bottom */ \ - - coeff4 * att * x[grid.index(i_r, i_theta + 1)] /* Top */ \ - ); \ - } \ - else { \ - /* i_theta % 2 == 0 and i_r % 2 == 0 */ \ - /* ---------- */ \ - /* o o o */ \ - /* ---------- */ \ - /* o X o */ \ - /* ---------- */ \ - /* o o o */ \ - /* ---------- */ \ - /* Fill temp(i-1,j) */ \ - temp[grid.index(i_r - 1, i_theta)] -= \ - (-coeff1 * arr * x[grid.index(i_r, i_theta)] /* Right */ \ - - 0.25 * art * x[grid.index(i_r, i_theta + 1)] /* Top Right */ \ - + 0.25 * art * x[grid.index(i_r, i_theta - 1)]); /* Bottom Right */ \ - /* Fill temp(i+1,j) */ \ - temp[grid.index(i_r + 1, i_theta)] -= \ - (-coeff2 * arr * x[grid.index(i_r, i_theta)] /* Left */ \ - + 0.25 * art * x[grid.index(i_r, i_theta + 1)] /* Top Left */ \ - - 0.25 * art * x[grid.index(i_r, i_theta - 1)]); /* Bottom Left */ \ - } \ - } \ - } \ - /* --------------------- */ \ - /* Outside Section Parts */ \ - else if (node_color != smoother_color) { \ - if (i_theta & 1) { \ - if (i_r & 1) { \ - /* i_theta % 2 == 1 and i_r % 2 == 1 */ \ - /* ---------- */ \ - /* x o x */ \ - /* ---------- */ \ - /* o O o */ \ - /* ---------- */ \ - /* x o x */ \ - /* ---------- */ \ - /* Fill temp(i,j-1) */ \ - temp[grid.index(i_r, i_theta - 1)] -= \ - (-coeff3 * att * x[grid.index(i_r, i_theta)] /* Top */ \ - - 0.25 * art * x[grid.index(i_r + 1, i_theta)] /* Top Right */ \ - + 0.25 * art * x[grid.index(i_r - 1, i_theta)]); /* Top Left */ \ - /* Fill temp(i,j+1) */ \ - temp[grid.index(i_r, i_theta + 1)] -= \ - (-coeff4 * att * x[grid.index(i_r, i_theta)] /* Bottom */ \ - + 0.25 * art * x[grid.index(i_r + 1, i_theta)] /* Bottom Right */ \ - - 0.25 * art * x[grid.index(i_r - 1, i_theta)]); /* Bottom Left */ \ - } \ - else { \ - /* i_theta % 2 == 1 and i_r % 2 == 0 */ \ - /* ---------- */ \ - /* o x o */ \ - /* ---------- */ \ - /* o O o */ \ - /* ---------- */ \ - /* o x o */ \ - /* ---------- */ \ - \ - /* Nothing to do! */ \ - } \ - } \ - else { \ - if (i_r & 1) { \ - /* i_theta % 2 == 0 and i_r % 2 == 1 */ \ - /* ---------- */ \ - /* o o o */ \ - /* ---------- */ \ - /* x O x */ \ - /* ---------- */ \ - /* o o o */ \ - /* ---------- */ \ - /* Fill temp(i,j-1) */ \ - temp[grid.index(i_r, i_theta - 1)] -= \ - (-coeff3 * att * x[grid.index(i_r, i_theta)] /* Top */ \ - - 0.25 * art * x[grid.index(i_r + 1, i_theta)] /* Top Right */ \ - + 0.25 * art * x[grid.index(i_r - 1, i_theta)]); /* Top Left */ \ - /* Fill temp(i,j+1) */ \ - temp[grid.index(i_r, i_theta + 1)] -= \ - (-coeff4 * att * x[grid.index(i_r, i_theta)] /* Bottom */ \ - + 0.25 * art * x[grid.index(i_r + 1, i_theta)] /* Bottom Right */ \ - - 0.25 * art * x[grid.index(i_r - 1, i_theta)]); /* Bottom Left */ \ - } \ - else { \ - /* i_theta % 2 == 0 and i_r % 2 == 0 */ \ - /* ---------- */ \ - /* o o o */ \ - /* ---------- */ \ - /* o X o */ \ - /* ---------- */ \ - /* o o o */ \ - /* ---------- */ \ - /* Fill temp(i,j-1) */ \ - temp[grid.index(i_r, i_theta - 1)] -= \ - (-coeff3 * att * x[grid.index(i_r, i_theta)] /* Top */ \ - - 0.25 * art * x[grid.index(i_r + 1, i_theta)] /* Top Right */ \ - + 0.25 * art * x[grid.index(i_r - 1, i_theta)]); /* Top Left */ \ - /* Fill temp(i,j+1) */ \ - temp[grid.index(i_r, i_theta + 1)] -= \ - (-coeff4 * att * x[grid.index(i_r, i_theta)] /* Bottom */ \ - + 0.25 * art * x[grid.index(i_r + 1, i_theta)] /* Bottom Right */ \ - - 0.25 * art * x[grid.index(i_r - 1, i_theta)]); /* Bottom Left */ \ - } \ - } \ - } \ - } \ - else if (i_r == grid.numberSmootherCircles() - 1) { \ - double h1 = grid.radialSpacing(i_r - 1); \ - double h2 = grid.radialSpacing(i_r); \ - double k1 = grid.angularSpacing(i_theta - 1); \ - double k2 = grid.angularSpacing(i_theta); \ - double coeff1 = 0.5 * (k1 + k2) / h1; \ - double coeff2 = 0.5 * (k1 + k2) / h2; \ - double coeff3 = 0.5 * (h1 + h2) / k1; \ - double coeff4 = 0.5 * (h1 + h2) / k2; \ - /* -------------------- */ \ - /* Inside Section Parts */ \ - if (node_color == smoother_color) { \ - /* Dont give to the right when this case occurs! */ \ - /* i_theta % 2 = 0 and i_r % 2 == 1 */ \ - /* | o | o | o || o o o o */ \ - /* | | | || -------------- */ \ - /* | o | x | O || x o x o */ \ - /* | | | || -------------- */ \ - /* | o | o | o || o o o o */ \ - if ((!(i_r & 1) || (i_theta & 1))) { \ - /* Fill temp(i+1,j) */ \ - temp[grid.index(i_r + 1, i_theta)] -= \ - (-coeff2 * arr * x[grid.index(i_r, i_theta)] /* Left */ \ - + 0.25 * art * x[grid.index(i_r, i_theta + 1)] /* Top Left */ \ - - 0.25 * art * x[grid.index(i_r, i_theta - 1)]); /* Bottom Left */ \ - } \ - } \ - /* --------------------- */ \ - /* Outside Section Parts */ \ - else if (node_color != smoother_color) { \ - /* Nothing to be done here */ \ - } \ - } \ - else if (i_r == grid.numberSmootherCircles()) { \ - double h1 = grid.radialSpacing(i_r - 1); \ - double h2 = grid.radialSpacing(i_r); \ - double k1 = grid.angularSpacing(i_theta - 1); \ - double k2 = grid.angularSpacing(i_theta); \ - double coeff1 = 0.5 * (k1 + k2) / h1; \ - double coeff2 = 0.5 * (k1 + k2) / h2; \ - double coeff3 = 0.5 * (h1 + h2) / k1; \ - double coeff4 = 0.5 * (h1 + h2) / k2; \ - /* -------------------- */ \ - /* Inside Section Parts */ \ - if (node_color == smoother_color) { \ - if (i_theta & 1) { \ - if (i_r & 1) { \ - /* i_theta % 2 == 1 and i_r % 2 == 1 */ \ - /* | x | o | x || o x o x */ \ - /* | | | || -------------- */ \ - /* | o | o | o || O o o o */ \ - /* | | | || -------------- */ \ - /* | x | o | x || o x o x */ \ - /* Fill temp(i,j) */ \ - temp[grid.index(i_r, i_theta)] -= \ - (-coeff1 * arr * x[grid.index(i_r - 1, i_theta)] /* Left */ \ - - coeff3 * att * x[grid.index(i_r, i_theta - 1)] /* Bottom */ \ - - coeff4 * att * x[grid.index(i_r, i_theta + 1)] /* Top */ \ - ); \ - /* Fill temp(i+1,j) */ \ - temp[grid.index(i_r + 1, i_theta)] -= \ - (+0.25 * art * x[grid.index(i_r, i_theta + 1)] /* Top Left */ \ - - 0.25 * art * x[grid.index(i_r, i_theta - 1)]); /* Bottom Left */ \ - } \ - else { \ - /* i_theta % 2 == 1 and i_r % 2 == 0 */ \ - /* | o | x | o || x o x o */ \ - /* | | | || -------------- */ \ - /* | o | o | o || O o o o */ \ - /* | | | || -------------- */ \ - /* | o | x | o || x o x o */ \ - /* Fill temp(i,j) */ \ - temp[grid.index(i_r, i_theta)] -= \ - (-coeff1 * arr * x[grid.index(i_r - 1, i_theta)] /* Left */ \ - - coeff3 * att * x[grid.index(i_r, i_theta - 1)] /* Bottom */ \ - - coeff4 * att * x[grid.index(i_r, i_theta + 1)] /* Top */ \ - ); \ - /* Fill temp(i+1,j) */ \ - temp[grid.index(i_r + 1, i_theta)] -= \ - (+0.25 * art * x[grid.index(i_r, i_theta + 1)] /* Top Left */ \ - - 0.25 * art * x[grid.index(i_r, i_theta - 1)]); /* Bottom Left */ \ - } \ - } \ - else { \ - if (i_r & 1) { \ - /* i_theta % 2 == 0 and i_r % 2 == 1 */ \ - /* | o | o | o || o o o o */ \ - /* | | | || -------------- */ \ - /* | x | o | x || O x o x */ \ - /* | | | || -------------- */ \ - /* | o | o | o || o o o o */ \ - /* Fill temp(i,j) */ \ - temp[grid.index(i_r, i_theta)] -= \ - (-coeff1 * arr * x[grid.index(i_r - 1, i_theta)] /* Left */ \ - - coeff2 * arr * x[grid.index(i_r + 1, i_theta)] /* Right */ \ - - coeff3 * att * x[grid.index(i_r, i_theta - 1)] /* Bottom */ \ - - coeff4 * att * x[grid.index(i_r, i_theta + 1)] /* Top */ \ - ); \ - } \ - else { \ - /* i_theta % 2 == 0 and i_r % 2 == 0 */ \ - /* | o | o | o || o o o o */ \ - /* | | | || -------------- */ \ - /* | o | x | o || X o x o */ \ - /* | | | || -------------- */ \ - /* | o | o | o || o o o o */ \ - /* Fill temp(i+1,j) */ \ - temp[grid.index(i_r + 1, i_theta)] -= \ - (-coeff2 * arr * x[grid.index(i_r, i_theta)] /* Left */ \ - + 0.25 * art * x[grid.index(i_r, i_theta + 1)] /* Top Left */ \ - - 0.25 * art * x[grid.index(i_r, i_theta - 1)]); /* Bottom Left */ \ - } \ - } \ - } \ - /* --------------------- */ \ - /* Outside Section Parts */ \ - else if (node_color != smoother_color) { \ - /* Dont give to bottom and up when this case occurs! */ \ - /* i_theta % 2 == 1 and i_r % 2 == 0 */ \ - /* | o | x | o || x o x o */ \ - /* | | | || -------------- */ \ - /* | o | o | o || O o o o */ \ - /* | | | || -------------- */ \ - /* | o | x | o || x o x o */ \ - if (i_r & 1 || !(i_theta & 1)) { \ - /* Fill temp(i,j-1) */ \ - temp[grid.index(i_r, i_theta - 1)] -= \ - (-coeff3 * att * x[grid.index(i_r, i_theta)] /* Top */ \ - - 0.25 * art * x[grid.index(i_r + 1, i_theta)] /* Top Right */ \ - + 0.25 * art * x[grid.index(i_r - 1, i_theta)]); /* Top Left */ \ - /* Fill temp(i,j+1) */ \ - temp[grid.index(i_r, i_theta + 1)] -= \ - (-coeff4 * att * x[grid.index(i_r, i_theta)] /* Bottom */ \ - + 0.25 * art * x[grid.index(i_r + 1, i_theta)] /* Bottom Right */ \ - - 0.25 * art * x[grid.index(i_r - 1, i_theta)]); /* Bottom Left */ \ - } \ - } \ - } \ - else if (i_r == grid.nr() - 2) { \ - assert(i_r & 1); \ - \ - double h1 = grid.radialSpacing(i_r - 1); \ - double h2 = grid.radialSpacing(i_r); \ - double k1 = grid.angularSpacing(i_theta - 1); \ - double k2 = grid.angularSpacing(i_theta); \ - double coeff1 = 0.5 * (k1 + k2) / h1; \ - double coeff2 = 0.5 * (k1 + k2) / h2; \ - double coeff3 = 0.5 * (h1 + h2) / k1; \ - double coeff4 = 0.5 * (h1 + h2) / k2; \ - /* -------------------- */ \ - /* Inside Section Parts */ \ - if (node_color == smoother_color) { \ - if (i_theta & 1) { \ - /* i_theta % 2 == 1 */ \ - /* ---------------|| */ \ - /* o x o x || */ \ - /* ---------------|| */ \ - /* o o O o || */ \ - /* ---------------|| */ \ - /* o x o x || */ \ - /* ---------------|| */ \ - /* Fill temp(i,j) */ \ - temp[grid.index(i_r, i_theta)] -= (-coeff3 * att * x[grid.index(i_r, i_theta - 1)] /* Bottom */ \ - - coeff4 * att * x[grid.index(i_r, i_theta + 1)] /* Top */ \ - ); \ - /* Fill temp(i-1,j) */ \ - temp[grid.index(i_r - 1, i_theta)] -= \ - (-0.25 * art * x[grid.index(i_r, i_theta + 1)] /* Top Right */ \ - + 0.25 * art * x[grid.index(i_r, i_theta - 1)]); /* Bottom Right */ \ - \ - /* "Right" is part of the radial Asc smoother matrices, */ \ - /* but is shifted over to the rhs to make the radial Asc smoother matrices symmetric. */ \ - /* Note that the circle Asc smoother matrices are symmetric by default. */ \ - /* Note that rhs[grid.index(i_r + 1, i_theta)] contains the correct boundary value of u_D. */ \ - temp[grid.index(i_r, i_theta)] -= /* Right: Symmetry shift! */ \ - -coeff2 * arr * rhs[grid.index(i_r + 1, i_theta)]; \ - } \ - else { \ - /* ---------------|| */ \ - /* o o o o || */ \ - /* ---------------|| */ \ - /* o x O x || */ \ - /* ---------------|| */ \ - /* o o o o || */ \ - /* ---------------|| */ \ - /* Fill temp(i,j) */ \ - temp[grid.index(i_r, i_theta)] -= (-coeff1 * arr * x[grid.index(i_r - 1, i_theta)] /* Left */ \ - - coeff2 * arr * x[grid.index(i_r + 1, i_theta)] /* Right */ \ - - coeff3 * att * x[grid.index(i_r, i_theta - 1)] /* Bottom */ \ - - coeff4 * att * x[grid.index(i_r, i_theta + 1)] /* Top */ \ - ); \ - } \ - } \ - /* --------------------- */ \ - /* Outside Section Parts */ \ - else if (node_color != smoother_color) { \ - /* Fill temp(i,j-1) */ \ - temp[grid.index(i_r, i_theta - 1)] -= (-coeff3 * att * x[grid.index(i_r, i_theta)] /* Top */ \ - - 0.25 * art * x[grid.index(i_r + 1, i_theta)] /* Top Right */ \ - + 0.25 * art * x[grid.index(i_r - 1, i_theta)]); /* Top Left */ \ - /* Fill temp(i,j+1) */ \ - temp[grid.index(i_r, i_theta + 1)] -= \ - (-coeff4 * att * x[grid.index(i_r, i_theta)] /* Bottom */ \ - + 0.25 * art * x[grid.index(i_r + 1, i_theta)] /* Bottom Right */ \ - - 0.25 * art * x[grid.index(i_r - 1, i_theta)]); /* Bottom Left */ \ - } \ - } \ - else if (i_r == grid.nr() - 1) { \ - assert(!(i_r & 1)); \ - \ - double h1 = grid.radialSpacing(i_r - 1); \ - double k1 = grid.angularSpacing(i_theta - 1); \ - double k2 = grid.angularSpacing(i_theta); \ - double coeff1 = 0.5 * (k1 + k2) / h1; \ - /* -------------------- */ \ - /* Inside Section Parts */ \ - if (node_color == smoother_color) { \ - if (i_theta & 1) { \ - /* i_theta % 2 == 1 */ \ - /* -----------|| */ \ - /* x o x || */ \ - /* -----------|| */ \ - /* o o O || */ \ - /* -----------|| */ \ - /* x o x || */ \ - /* -----------|| */ \ - /* Fill temp(i-1,j) */ \ - temp[grid.index(i_r - 1, i_theta)] -= \ - (-0.25 * art * x[grid.index(i_r, i_theta + 1)] /* Top Right */ \ - + 0.25 * art * x[grid.index(i_r, i_theta - 1)] /* Bottom Right */ \ - ); \ - /* "Right" is part of the radial Asc smoother matrices, */ \ - /* but is shifted over to the rhs to make the radial Asc smoother matrices symmetric. */ \ - /* Note that the circle Asc smoother matrices are symmetric by default. */ \ - /* Note that rhs[grid.index(i_r, i_theta)] contains the correct boundary value of u_D. */ \ - temp[grid.index(i_r - 1, i_theta)] -= (-coeff1 * arr * rhs[grid.index(i_r, i_theta)] /* Right */ \ - ); \ - } \ - else { \ - /* -----------|| */ \ - /* o o o || */ \ - /* -----------|| */ \ - /* x o X || */ \ - /* -----------|| */ \ - /* o o o || */ \ - /* -----------|| */ \ - /* Fill temp(i-1,j) */ \ - temp[grid.index(i_r - 1, i_theta)] -= \ - (-coeff1 * arr * x[grid.index(i_r, i_theta)] /* Right */ \ - - 0.25 * art * x[grid.index(i_r, i_theta + 1)] /* Top Right */ \ - + 0.25 * art * x[grid.index(i_r, i_theta - 1)]); /* Bottom Right */ \ - } \ - } \ - /* --------------------- */ \ - /* Outside Section Parts */ \ - else if (node_color != smoother_color) { \ - /* Nothing to be done here */ \ - } \ - } \ - } while (0) +inline void nodeApplyAscOrthoCircleGive(int i_r, int i_theta, const PolarGrid& grid, bool DirBC_Interior, + SmootherColor smoother_color, ConstVector& x, ConstVector& rhs, + Vector& result, double arr, double att, double art, double detDF, + double coeff_beta) +{ + assert(i_r >= 0 && i_r <= grid.numberSmootherCircles()); + + bool isOddNumberSmootherCircles = (grid.numberSmootherCircles() & 1); + bool isOddRadialIndex = (i_r & 1); + + SmootherColor node_color = + (isOddNumberSmootherCircles == isOddRadialIndex) ? SmootherColor::White : SmootherColor::Black; + + /* -------------------- */ + /* Node in the interior */ + /* -------------------- */ + if (i_r > 0 && i_r < grid.numberSmootherCircles()) { + double h1 = grid.radialSpacing(i_r - 1); + double h2 = grid.radialSpacing(i_r); + double k1 = grid.angularSpacing(i_theta - 1); + double k2 = grid.angularSpacing(i_theta); + + double coeff1 = 0.5 * (k1 + k2) / h1; + double coeff2 = 0.5 * (k1 + k2) / h2; + double coeff3 = 0.5 * (h1 + h2) / k1; + double coeff4 = 0.5 * (h1 + h2) / k2; + + int i_theta_M1 = grid.wrapThetaIndex(i_theta - 1); + int i_theta_P1 = grid.wrapThetaIndex(i_theta + 1); + + int center = grid.index(i_r, i_theta); + int left = grid.index(i_r - 1, i_theta); + int right = grid.index(i_r + 1, i_theta); + int bottom = grid.index(i_r, i_theta_M1); + int top = grid.index(i_r, i_theta_P1); + + /* -------------------- */ + /* Inside Section Parts */ + if (node_color == smoother_color) { + if (i_r & 1) { + if (i_theta & 1) { + /* i_r % 2 == 1 and i_theta % 2 == 1 */ + /* | x | o | x | */ + /* | | | | */ + /* | o | O | o | */ + /* | | | | */ + /* | x | o | x | */ + + /* Fill result(i,j) */ + result[center] -= (-coeff1 * arr * x[left] /* Left */ + - coeff2 * arr * x[right] /* Right */ + ); /* Fill result(i,j-1) */ + result[bottom] -= (-0.25 * art * x[right] /* Top Right */ + + 0.25 * art * x[left]); /* Top Left */ + /* Fill result(i,j+1) */ + result[top] -= (+0.25 * art * x[right] /* Bottom Right */ + - 0.25 * art * x[left]); /* Bottom Left */ + } + else { + /* i_r % 2 == 1 and i_theta % 2 == 0 */ + /* | o | o | o | */ + /* | | | | */ + /* | x | O | x | */ + /* | | | | */ + /* | o | o | o | */ + + /* Fill result(i,j) */ + result[center] -= (-coeff1 * arr * x[left] /* Left */ + - coeff2 * arr * x[right] /* Right */ + ); + /* Fill result(i,j-1) */ + result[bottom] -= (-0.25 * art * x[right] /* Top Right */ + + 0.25 * art * x[left]); /* Top Left */ + /* Fill result(i,j+1) */ + result[top] -= (+0.25 * art * x[right] /* Bottom Right */ + - 0.25 * art * x[left]); /* Bottom Left */ + } + } + else { + if (i_theta & 1) { + /* i_r % 2 == 0 and i_theta % 2 == 1 */ + /* | o | x | o | */ + /* | | | | */ + /* | o | O | o | */ + /* | | | | */ + /* | o | x | o | */ + + /* Fill result(i,j) */ + result[center] -= (-coeff1 * arr * x[left] /* Left */ + - coeff2 * arr * x[right] /* Right */ + - coeff3 * att * x[bottom] /* Bottom */ + - coeff4 * att * x[top] /* Top */ + ); + } + else { + /* i_r % 2 == 0 and i_theta % 2 == 0 */ + /* | o | o | o | */ + /* | | | | */ + /* | o | X | o | */ + /* | | | | */ + /* | o | o | o | */ + + /* Fill result(i,j-1) */ + result[bottom] -= (-coeff3 * att * x[center] /* Top */ + - 0.25 * art * x[right] /* Top Right */ + + 0.25 * art * x[left]); /* Top Left */ + /* Fill result(i,j+1) */ + result[top] -= (-coeff4 * att * x[center] /* Bottom */ + + 0.25 * art * x[right] /* Bottom Right */ + - 0.25 * art * x[left]); /* Bottom Left */ + } + } + } + /* --------------------- */ + /* Outside Section Parts */ + else if (node_color != smoother_color) { + if (i_r & 1) { + if (i_theta & 1) { + /* i_r % 2 == 1 and i_theta % 2 == 1 */ + /* | x | o | x | */ + /* | | | | */ + /* | o | O | o | */ + /* | | | | */ + /* | x | o | x | */ + + /* Fill result(i-1,j) */ + if (i_r > 1 || !DirBC_Interior) { + result[left] -= (-coeff1 * arr * x[center] /* Right */ + - 0.25 * art * x[top] /* Top Right */ + + 0.25 * art * x[bottom]); /* Bottom Right */ + } /* Fill result(i+1,j) */ + if (i_r < grid.numberSmootherCircles() - 1) { + result[right] -= (-coeff2 * arr * x[center] /* Left */ + + 0.25 * art * x[top] /* Top Left */ + - 0.25 * art * x[bottom]); /* Bottom Left */ + } + } + else { + /* i_r % 2 == 1 and i_theta % 2 == 0 */ + /* | o | o | o | */ + /* | | | | */ + /* | x | O | x | */ + /* | | | | */ + /* | o | o | o | */ + + /* Nothing to do! */ + } + } + else { + if (i_theta & 1) { + /* i_r % 2 == 0 and i_theta % 2 == 1 */ + /* | o | x | o | */ + /* | | | | */ + /* | o | O | o | */ + /* | | | | */ + /* | o | x | o | */ + + /* Fill result(i-1,j) */ + if (i_r > 1 || !DirBC_Interior) { + result[left] -= (-coeff1 * arr * x[center] /* Right */ + - 0.25 * art * x[top] /* Top Right */ + + 0.25 * art * x[bottom]); /* Bottom Right */ + } + /* Fill result(i+1,j) */ + if (i_r < grid.numberSmootherCircles() - 1) { + result[right] -= (-coeff2 * arr * x[center] /* Left */ + + 0.25 * art * x[top] /* Top Left */ + - 0.25 * art * x[bottom]); /* Bottom Left */ + } + } + else { + /* i_r % 2 == 0 and i_theta % 2 == 0 */ + /* | o | o | o | */ + /* | | | | */ + /* | o | X | o | */ + /* | | | | */ + /* | o | o | o | */ + + /* Fill result(i-1,j) */ + if (i_r > 1 || !DirBC_Interior) { + result[left] -= (-coeff1 * arr * x[center] /* Right */ + - 0.25 * art * x[top] /* Top Right */ + + 0.25 * art * x[bottom]); /* Bottom Right */ + } /* Fill result(i+1,j) */ + if (i_r < grid.numberSmootherCircles() - 1) { + result[right] -= (-coeff2 * arr * x[center] /* Left */ + + 0.25 * art * x[top] /* Top Left */ + - 0.25 * art * x[bottom]); /* Bottom Left */ + } + } + } + } + } + /* -------------------- */ + /* Node on the boundary */ + /* -------------------- */ + else if (i_r == 0) { + /* ------------------------------------------------ */ + /* Case 1: Dirichlet boundary on the inner boundary */ + /* ------------------------------------------------ */ + if (DirBC_Interior) { + double h2 = grid.radialSpacing(i_r); + double k1 = grid.angularSpacing(i_theta - 1); + double k2 = grid.angularSpacing(i_theta); + + double coeff2 = 0.5 * (k1 + k2) / h2; + + int i_theta_M1 = grid.wrapThetaIndex(i_theta - 1); + int i_theta_P1 = grid.wrapThetaIndex(i_theta + 1); + + int center = grid.index(i_r, i_theta); + int right = grid.index(i_r + 1, i_theta); + int bottom = grid.index(i_r, i_theta_M1); + int top = grid.index(i_r, i_theta_P1); + + /* -------------------- */ + /* Inside Section Parts */ + if (node_color == smoother_color) { + if (i_theta & 1) { + /* i_theta % 2 == 1 */ + /* || x | o | x | */ + /* || | | | */ + /* || O | o | o | */ + /* || | | | */ + /* || x | o | x | */ + + /* Nothing to do! */ + } + else { + /* i_theta % 2 == 0 */ + /* || o | o | o | */ + /* || | | | */ + /* || X | o | x | */ + /* || | | | */ + /* || o | o | o | */ + + /* Nothing to do! */ + } + } + /* --------------------- */ + /* Outside Section Parts */ + else if (node_color != smoother_color) { + /* Fill result(i+1,j) */ + result[right] -= (-coeff2 * arr * x[center] /* Left */ + + 0.25 * art * x[top] /* Top Left */ + - 0.25 * art * x[bottom]); /* Bottom Left */ + } + } + else { + /* ------------------------------------------------------------- */ + /* Case 2: Across origin discretization on the interior boundary */ + /* ------------------------------------------------------------- */ + // h1 gets replaced with 2 * R0. + // (i_r-1,i_theta) gets replaced with (i_r, i_theta + (grid.ntheta()/2)). + // Some more adjustments from the changing the 9-point stencil to the artifical 7-point stencil. + double h1 = 2.0 * grid.radius(0); + double h2 = grid.radialSpacing(i_r); + double k1 = grid.angularSpacing(i_theta - 1); + double k2 = grid.angularSpacing(i_theta); + + double coeff1 = 0.5 * (k1 + k2) / h1; + double coeff2 = 0.5 * (k1 + k2) / h2; + double coeff3 = 0.5 * (h1 + h2) / k1; + double coeff4 = 0.5 * (h1 + h2) / k2; + + const int i_theta_M1 = grid.wrapThetaIndex(i_theta - 1); + const int i_theta_P1 = grid.wrapThetaIndex(i_theta + 1); + const int i_theta_Across = grid.wrapThetaIndex(i_theta + grid.ntheta() / 2); + + const int left = grid.index(i_r, i_theta_Across); + const int bottom = grid.index(i_r, i_theta_M1); + const int center = grid.index(i_r, i_theta); + const int top = grid.index(i_r, i_theta_P1); + const int right = grid.index(i_r + 1, i_theta); + + /* -------------------- */ + /* Inside Section Parts */ + if (node_color == smoother_color) { + if (i_theta & 1) { + /* i_theta % 2 == 1 */ + /* -| x | o | x | */ + /* -| | | | */ + /* -| O | o | o | */ + /* -| | | | */ + /* -| x | o | x | */ + + /* Fill result(i,j) */ + result[center] -= ( + /* - coeff1 * arr * x[left] // Left: Not in Asc_ortho */ + -coeff2 * arr * x[right] /* Right */ + - coeff3 * att * x[bottom] /* Bottom */ + - coeff4 * att * x[top] /* Top */ + ); + } + else { + /* i_theta % 2 == 0 */ + /* -| o | o | o | */ + /* -| | | | */ + /* -| X | o | x | */ + /* -| | | | */ + /* -| o | o | o | */ + + /* Fill result(i,j-1) */ + result[bottom] -= (-coeff3 * att * x[center] /* Top */ + - 0.25 * art * x[right]); /* Top Right */ + /* + 0.25 * art * x[grid.index(i_r-1,i_theta)]; // Top Left: REMOVED DUE TO ARTIFICAL 7 POINT STENCIL */ + /* Fill result(i,j+1) */ + result[top] -= (-coeff4 * att * x[center] /* Bottom */ + + 0.25 * art * x[right]); /* Bottom Right */ + /* - 0.25 * art * x[grid.index(i_r-1,i_theta)]; // Bottom Left: REMOVED DUE TO ARTIFICAL 7 POINT STENCIL */ + } + } + /* --------------------- */ + /* Outside Section Parts */ + else if (node_color != smoother_color) { + /* Fill result(i+1,j) */ + result[right] -= (-coeff2 * arr * x[center] /* Left */ + + 0.25 * art * x[top] /* Top Left */ + - 0.25 * art * x[bottom]); /* Bottom Left */ + } + } + } + /* ----------------------------- */ + /* Node next to circular section */ + /* ----------------------------- */ + else if (i_r == grid.numberSmootherCircles()) { + assert(node_color == SmootherColor::White); + + if (smoother_color == SmootherColor::Black) { + double h1 = grid.radialSpacing(i_r - 1); + double h2 = grid.radialSpacing(i_r); + double k1 = grid.angularSpacing(i_theta - 1); + double k2 = grid.angularSpacing(i_theta); + + double coeff1 = 0.5 * (k1 + k2) / h1; + double coeff2 = 0.5 * (k1 + k2) / h2; + double coeff3 = 0.5 * (h1 + h2) / k1; + double coeff4 = 0.5 * (h1 + h2) / k2; + + int i_theta_M1 = grid.wrapThetaIndex(i_theta - 1); + int i_theta_P1 = grid.wrapThetaIndex(i_theta + 1); + + int center = grid.index(i_r, i_theta); + int left = grid.index(i_r - 1, i_theta); + int bottom = grid.index(i_r, i_theta_M1); + int top = grid.index(i_r, i_theta_P1); + + /* i_theta % 2 == 1 and i_r % 2 == 1 */ + /* | x | o | x || o x o x */ + /* | | | || -------------- */ + /* | o | o | o || O o o o */ + /* | | | || -------------- */ + /* | x | o | x || o x o x */ + /* -> Give Left */ + + /* i_theta % 2 == 1 and i_r % 2 == 0 */ + /* | o | x | o || x o x o */ + /* | | | || -------------- */ + /* | o | o | o || O o o o */ + /* | | | || -------------- */ + /* | o | x | o || x o x o */ + /* -> Give Left */ + + /* i_theta % 2 == 0 and i_r % 2 == 1 */ + /* | o | o | o || o o o o */ + /* | | | || -------------- */ + /* | x | o | x || O x o x */ + /* | | | || -------------- */ + /* | o | o | o || o o o o */ + /* -> Don't give to the Left! */ + + /* i_theta % 2 == 0 and i_r % 2 == 0 */ + /* | o | o | o || o o o o */ + /* | | | || -------------- */ + /* | o | x | o || X o x o */ + /* | | | || -------------- */ + /* | o | o | o || o o o o */ + /* -> Give Left */ + + if (i_theta & 1 || !(i_r & 1)) { + /* --------------------- */ + /* Outside Section Parts */ + + /* Fill result(i-1,j) */ + result[left] -= (-coeff1 * arr * x[center] /* Right */ + - 0.25 * art * x[top] /* Top Right */ + + 0.25 * art * x[bottom]); /* Bottom Right */ + } + } + } +} + +inline void nodeApplyAscOrthoRadialGive(int i_r, int i_theta, const PolarGrid& grid, bool DirBC_Interior, + SmootherColor smoother_color, ConstVector& x, ConstVector& rhs, + Vector& result, double arr, double att, double art, double detDF, + double coeff_beta) +{ + assert(i_r >= grid.numberSmootherCircles() - 1 && i_r < grid.nr()); + + SmootherColor node_color = (i_theta & 1) ? SmootherColor::White : SmootherColor::Black; + + /* -------------------- */ + /* Node in the interior */ + /* -------------------- */ + if (i_r > grid.numberSmootherCircles() && i_r < grid.nr() - 2) { + double h1 = grid.radialSpacing(i_r - 1); + double h2 = grid.radialSpacing(i_r); + double k1 = grid.angularSpacing(i_theta - 1); + double k2 = grid.angularSpacing(i_theta); + + double coeff1 = 0.5 * (k1 + k2) / h1; + double coeff2 = 0.5 * (k1 + k2) / h2; + double coeff3 = 0.5 * (h1 + h2) / k1; + double coeff4 = 0.5 * (h1 + h2) / k2; + + int i_theta_M1 = grid.wrapThetaIndex(i_theta - 1); + int i_theta_P1 = grid.wrapThetaIndex(i_theta + 1); + + int center = grid.index(i_r, i_theta); + int left = grid.index(i_r - 1, i_theta); + int right = grid.index(i_r + 1, i_theta); + int bottom = grid.index(i_r, i_theta_M1); + int top = grid.index(i_r, i_theta_P1); + + /* -------------------- */ + /* Inside Section Parts */ + if (node_color == smoother_color) { + if (i_theta & 1) { + if (i_r & 1) { + /* i_theta % 2 == 1 and i_r % 2 == 1 */ + /* ---------- */ + /* x o x */ + /* ---------- */ + /* o O o */ + /* ---------- */ + /* x o x */ + /* ---------- */ + + /* Fill result(i,j) */ + result[center] -= (-coeff3 * att * x[bottom] /* Bottom */ + - coeff4 * att * x[top] /* Top */ + ); + /* Fill result(i-1,j) */ + result[left] -= (-0.25 * art * x[top] /* Top Right */ + + 0.25 * art * x[bottom]); /* Bottom Right */ + /* Fill result(i+1,j) */ + result[right] -= (+0.25 * art * x[top] /* Top Left */ + - 0.25 * art * x[bottom]); /* Bottom Left */ + } + else { + /* i_theta % 2 == 1 and i_r % 2 == 0 */ + /* ---------- */ + /* o x o */ + /* ---------- */ + /* o O o */ + /* ---------- */ + /* o x o */ + /* ---------- */ + + /* Fill result(i,j) */ + result[center] -= (-coeff3 * att * x[bottom] /* Bottom */ + - coeff4 * att * x[top] /* Top */ + ); + /* Fill result(i-1,j) */ + result[left] -= (-0.25 * art * x[top] /* Top Right */ + + 0.25 * art * x[bottom]); /* Bottom Right */ + /* Fill result(i+1,j) */ + result[right] -= (+0.25 * art * x[top] /* Top Left */ + - 0.25 * art * x[bottom]); /* Bottom Left */ + } + } + else { + if (i_r & 1) { + /* i_theta % 2 == 0 and i_r % 2 == 1 */ + /* ---------- */ + /* o o o */ + /* ---------- */ + /* x O x */ + /* ---------- */ + /* o o o */ + /* ---------- */ + + /* Fill result(i,j) */ + result[center] -= (-coeff1 * arr * x[left] /* Left */ + - coeff2 * arr * x[right] /* Right */ + - coeff3 * att * x[bottom] /* Bottom */ + - coeff4 * att * x[top] /* Top */ + ); + } + else { + /* i_theta % 2 == 0 and i_r % 2 == 0 */ + /* ---------- */ + /* o o o */ + /* ---------- */ + /* o X o */ + /* ---------- */ + /* o o o */ + /* ---------- */ + + /* Fill result(i-1,j) */ + result[left] -= (-coeff1 * arr * x[center] /* Right */ + - 0.25 * art * x[top] /* Top Right */ + + 0.25 * art * x[bottom]); /* Bottom Right */ + /* Fill result(i+1,j) */ + result[right] -= (-coeff2 * arr * x[center] /* Left */ + + 0.25 * art * x[top] /* Top Left */ + - 0.25 * art * x[bottom]); /* Bottom Left */ + } + } + } + /* --------------------- */ + /* Outside Section Parts */ + else if (node_color != smoother_color) { + if (i_theta & 1) { + if (i_r & 1) { + /* i_theta % 2 == 1 and i_r % 2 == 1 */ + /* ---------- */ + /* x o x */ + /* ---------- */ + /* o O o */ + /* ---------- */ + /* x o x */ + /* ---------- */ + + /* Fill result(i,j-1) */ + result[bottom] -= (-coeff3 * att * x[center] /* Top */ + - 0.25 * art * x[right] /* Top Right */ + + 0.25 * art * x[left]); /* Top Left */ + /* Fill result(i,j+1) */ + result[top] -= (-coeff4 * att * x[center] /* Bottom */ + + 0.25 * art * x[right] /* Bottom Right */ + - 0.25 * art * x[left]); /* Bottom Left */ + } + else { + /* i_theta % 2 == 1 and i_r % 2 == 0 */ + /* ---------- */ + /* o x o */ + /* ---------- */ + /* o O o */ + /* ---------- */ + /* o x o */ + /* ---------- */ + + /* Nothing to do! */ + } + } + else { + if (i_r & 1) { + /* i_theta % 2 == 0 and i_r % 2 == 1 */ + /* ---------- */ + /* o o o */ + /* ---------- */ + /* x O x */ + /* ---------- */ + /* o o o */ + /* ---------- */ + + /* Fill result(i,j-1) */ + result[bottom] -= (-coeff3 * att * x[center] /* Top */ + - 0.25 * art * x[right] /* Top Right */ + + 0.25 * art * x[left]); /* Top Left */ + /* Fill result(i,j+1) */ + result[top] -= (-coeff4 * att * x[center] /* Bottom */ + + 0.25 * art * x[right] /* Bottom Right */ + - 0.25 * art * x[left]); /* Bottom Left */ + } + else { + /* i_theta % 2 == 0 and i_r % 2 == 0 */ + /* ---------- */ + /* o o o */ + /* ---------- */ + /* o X o */ + /* ---------- */ + /* o o o */ + /* ---------- */ + + /* Fill result(i,j-1) */ + result[bottom] -= (-coeff3 * att * x[center] /* Top */ + - 0.25 * art * x[right] /* Top Right */ + + 0.25 * art * x[left]); /* Top Left */ + /* Fill result(i,j+1) */ + result[top] -= (-coeff4 * att * x[center] /* Bottom */ + + 0.25 * art * x[right] /* Bottom Right */ + - 0.25 * art * x[left]); /* Bottom Left */ + } + } + } + } + else if (i_r == grid.numberSmootherCircles() - 1) { + double h1 = grid.radialSpacing(i_r - 1); + double h2 = grid.radialSpacing(i_r); + double k1 = grid.angularSpacing(i_theta - 1); + double k2 = grid.angularSpacing(i_theta); + + double coeff2 = 0.5 * (k1 + k2) / h2; + + int i_theta_M1 = grid.wrapThetaIndex(i_theta - 1); + int i_theta_P1 = grid.wrapThetaIndex(i_theta + 1); + + int center = grid.index(i_r, i_theta); + int right = grid.index(i_r + 1, i_theta); + int bottom = grid.index(i_r, i_theta_M1); + int top = grid.index(i_r, i_theta_P1); + + /* -------------------- */ + /* Inside Section Parts */ + if (node_color == smoother_color) { + /* Dont give to the right when this case occurs! */ + /* i_theta % 2 = 0 and i_r % 2 == 1 */ + /* | o | o | o || o o o o */ + /* | | | || -------------- */ + /* | o | x | O || x o x o */ + /* | | | || -------------- */ + /* | o | o | o || o o o o */ + if ((!(i_r & 1) || (i_theta & 1))) { + /* Fill result(i+1,j) */ + result[right] -= (-coeff2 * arr * x[center] /* Left */ + + 0.25 * art * x[top] /* Top Left */ + - 0.25 * art * x[bottom]); /* Bottom Left */ + } + } + /* --------------------- */ + /* Outside Section Parts */ + else if (node_color != smoother_color) { + /* Nothing to be done here */ + } + } + else if (i_r == grid.numberSmootherCircles()) { + double h1 = grid.radialSpacing(i_r - 1); + double h2 = grid.radialSpacing(i_r); + double k1 = grid.angularSpacing(i_theta - 1); + double k2 = grid.angularSpacing(i_theta); + + double coeff1 = 0.5 * (k1 + k2) / h1; + double coeff2 = 0.5 * (k1 + k2) / h2; + double coeff3 = 0.5 * (h1 + h2) / k1; + double coeff4 = 0.5 * (h1 + h2) / k2; + + int i_theta_M1 = grid.wrapThetaIndex(i_theta - 1); + int i_theta_P1 = grid.wrapThetaIndex(i_theta + 1); + + int center = grid.index(i_r, i_theta); + int left = grid.index(i_r - 1, i_theta); + int right = grid.index(i_r + 1, i_theta); + int bottom = grid.index(i_r, i_theta_M1); + int top = grid.index(i_r, i_theta_P1); + + /* -------------------- */ + /* Inside Section Parts */ + if (node_color == smoother_color) { + if (i_theta & 1) { + if (i_r & 1) { + /* i_theta % 2 == 1 and i_r % 2 == 1 */ + /* | x | o | x || o x o x */ + /* | | | || -------------- */ + /* | o | o | o || O o o o */ + /* | | | || -------------- */ + /* | x | o | x || o x o x */ + + /* Fill result(i,j) */ + result[center] -= (-coeff1 * arr * x[left] /* Left */ + - coeff3 * att * x[bottom] /* Bottom */ + - coeff4 * att * x[top] /* Top */ + ); + /* Fill result(i+1,j) */ + result[right] -= (+0.25 * art * x[top] /* Top Left */ + - 0.25 * art * x[bottom]); /* Bottom Left */ + } + else { + /* i_theta % 2 == 1 and i_r % 2 == 0 */ + /* | o | x | o || x o x o */ + /* | | | || -------------- */ + /* | o | o | o || O o o o */ + /* | | | || -------------- */ + /* | o | x | o || x o x o */ + + /* Fill result(i,j) */ + result[center] -= (-coeff1 * arr * x[left] /* Left */ + - coeff3 * att * x[bottom] /* Bottom */ + - coeff4 * att * x[top] /* Top */ + ); + /* Fill result(i+1,j) */ + result[right] -= (+0.25 * art * x[top] /* Top Left */ + - 0.25 * art * x[bottom]); /* Bottom Left */ + } + } + else { + if (i_r & 1) { + /* i_theta % 2 == 0 and i_r % 2 == 1 */ + /* | o | o | o || o o o o */ + /* | | | || -------------- */ + /* | x | o | x || O x o x */ + /* | | | || -------------- */ + /* | o | o | o || o o o o */ + + /* Fill result(i,j) */ + result[center] -= (-coeff1 * arr * x[left] /* Left */ + - coeff2 * arr * x[right] /* Right */ + - coeff3 * att * x[bottom] /* Bottom */ + - coeff4 * att * x[top] /* Top */ + ); + } + else { + /* i_theta % 2 == 0 and i_r % 2 == 0 */ + /* | o | o | o || o o o o */ + /* | | | || -------------- */ + /* | o | x | o || X o x o */ + /* | | | || -------------- */ + /* | o | o | o || o o o o */ + + /* Fill result(i+1,j) */ + result[right] -= (-coeff2 * arr * x[center] /* Left */ + + 0.25 * art * x[top] /* Top Left */ + - 0.25 * art * x[bottom]); /* Bottom Left */ + } + } + } + /* --------------------- */ + /* Outside Section Parts */ + else if (node_color != smoother_color) { + /* Dont give to bottom and up when this case occurs! */ + /* i_theta % 2 == 1 and i_r % 2 == 0 */ + /* | o | x | o || x o x o */ + /* | | | || -------------- */ + /* | o | o | o || O o o o */ + /* | | | || -------------- */ + /* | o | x | o || x o x o */ + if (i_r & 1 || !(i_theta & 1)) { + /* Fill result(i,j-1) */ + result[bottom] -= (-coeff3 * att * x[center] /* Top */ + - 0.25 * art * x[right] /* Top Right */ + + 0.25 * art * x[left]); /* Top Left */ + /* Fill result(i,j+1) */ + result[top] -= (-coeff4 * att * x[center] /* Bottom */ + + 0.25 * art * x[right] /* Bottom Right */ + - 0.25 * art * x[left]); /* Bottom Left */ + } + } + } + else if (i_r == grid.nr() - 2) { + assert(i_r & 1); + + double h1 = grid.radialSpacing(i_r - 1); + double h2 = grid.radialSpacing(i_r); + double k1 = grid.angularSpacing(i_theta - 1); + double k2 = grid.angularSpacing(i_theta); + + double coeff1 = 0.5 * (k1 + k2) / h1; + double coeff2 = 0.5 * (k1 + k2) / h2; + double coeff3 = 0.5 * (h1 + h2) / k1; + double coeff4 = 0.5 * (h1 + h2) / k2; + + int i_theta_M1 = grid.wrapThetaIndex(i_theta - 1); + int i_theta_P1 = grid.wrapThetaIndex(i_theta + 1); + + int center = grid.index(i_r, i_theta); + int left = grid.index(i_r - 1, i_theta); + int right = grid.index(i_r + 1, i_theta); + int bottom = grid.index(i_r, i_theta_M1); + int top = grid.index(i_r, i_theta_P1); + + /* -------------------- */ + /* Inside Section Parts */ + if (node_color == smoother_color) { + if (i_theta & 1) { + /* i_theta % 2 == 1 */ + /* ---------------|| */ + /* o x o x || */ + /* ---------------|| */ + /* o o O o || */ + /* ---------------|| */ + /* o x o x || */ + /* ---------------|| */ + + /* Fill result(i,j) */ + result[center] -= (-coeff3 * att * x[bottom] /* Bottom */ + - coeff4 * att * x[top] /* Top */ + ); + /* Fill result(i-1,j) */ + result[left] -= (-0.25 * art * x[top] /* Top Right */ + + 0.25 * art * x[bottom]); /* Bottom Right */ + + /* "Right" is part of the radial Asc smoother matrices, */ + /* but is shifted over to the rhs to make the radial Asc smoother matrices symmetric. */ + /* Note that the circle Asc smoother matrices are symmetric by default. */ + /* Note that rhs[right] contains the correct boundary value of u_D. */ + result[center] -= /* Right: Symmetry shift! */ + -coeff2 * arr * rhs[right]; + } + else { + /* ---------------|| */ + /* o o o o || */ + /* ---------------|| */ + /* o x O x || */ + /* ---------------|| */ + /* o o o o || */ + /* ---------------|| */ + + /* Fill result(i,j) */ + result[center] -= (-coeff1 * arr * x[left] /* Left */ + - coeff2 * arr * x[right] /* Right */ + - coeff3 * att * x[bottom] /* Bottom */ + - coeff4 * att * x[top] /* Top */ + ); + } + } + /* --------------------- */ + /* Outside Section Parts */ + else if (node_color != smoother_color) { + /* Fill result(i,j-1) */ + result[bottom] -= (-coeff3 * att * x[center] /* Top */ + - 0.25 * art * x[right] /* Top Right */ + + 0.25 * art * x[left]); /* Top Left */ + /* Fill result(i,j+1) */ + result[top] -= (-coeff4 * att * x[center] /* Bottom */ + + 0.25 * art * x[right] /* Bottom Right */ + - 0.25 * art * x[left]); /* Bottom Left */ + } + } + else if (i_r == grid.nr() - 1) { + assert(!(i_r & 1)); + + double h1 = grid.radialSpacing(i_r - 1); + double k1 = grid.angularSpacing(i_theta - 1); + double k2 = grid.angularSpacing(i_theta); + + double coeff1 = 0.5 * (k1 + k2) / h1; + + int i_theta_M1 = grid.wrapThetaIndex(i_theta - 1); + int i_theta_P1 = grid.wrapThetaIndex(i_theta + 1); + + int center = grid.index(i_r, i_theta); + int left = grid.index(i_r - 1, i_theta); + int bottom = grid.index(i_r, i_theta_M1); + int top = grid.index(i_r, i_theta_P1); + + /* -------------------- */ + /* Inside Section Parts */ + if (node_color == smoother_color) { + if (i_theta & 1) { + /* i_theta % 2 == 1 */ + /* -----------|| */ + /* x o x || */ + /* -----------|| */ + /* o o O || */ + /* -----------|| */ + /* x o x || */ + /* -----------|| */ + + /* Fill result(i-1,j) */ + result[left] -= (-0.25 * art * x[top] /* Top Right */ + + 0.25 * art * x[bottom] /* Bottom Right */ + ); + /* "Right" is part of the radial Asc smoother matrices, */ + /* but is shifted over to the rhs to make the radial Asc smoother matrices symmetric. */ + /* Note that the circle Asc smoother matrices are symmetric by default. */ + /* Note that rhs[center] contains the correct boundary value of u_D. */ + result[left] -= (-coeff1 * arr * rhs[center] /* Right */ + ); + } + else { + /* -----------|| */ + /* o o o || */ + /* -----------|| */ + /* x o X || */ + /* -----------|| */ + /* o o o || */ + /* -----------|| */ + + /* Fill result(i-1,j) */ + result[left] -= (-coeff1 * arr * x[center] /* Right */ + - 0.25 * art * x[top] /* Top Right */ + + 0.25 * art * x[bottom]); /* Bottom Right */ + } + } + /* --------------------- */ + /* Outside Section Parts */ + else if (node_color != smoother_color) { + /* Nothing to be done here */ + } + } +} void ExtrapolatedSmootherGive::applyAscOrthoCircleSection(const int i_r, const SmootherColor smoother_color, ConstVector x, ConstVector rhs, @@ -830,8 +901,8 @@ void ExtrapolatedSmootherGive::applyAscOrthoCircleSection(const int i_r, const S level_cache_.obtainValues(i_r, i_theta, index, r, theta, coeff_beta, arr, att, art, detDF); // Apply Asc Ortho at the current node - NODE_APPLY_ASC_ORTHO_CIRCLE_GIVE(i_r, i_theta, r, theta, grid_, DirBC_Interior_, smoother_color, x, rhs, temp, - arr, att, art, detDF, coeff_beta); + nodeApplyAscOrthoCircleGive(i_r, i_theta, grid_, DirBC_Interior_, smoother_color, x, rhs, temp, arr, att, art, + detDF, coeff_beta); } } @@ -850,8 +921,8 @@ void ExtrapolatedSmootherGive::applyAscOrthoRadialSection(const int i_theta, con level_cache_.obtainValues(i_r, i_theta, index, r, theta, coeff_beta, arr, att, art, detDF); // Apply Asc Ortho at the current node - NODE_APPLY_ASC_ORTHO_RADIAL_GIVE(i_r, i_theta, r, theta, grid_, DirBC_Interior_, smoother_color, x, rhs, temp, - arr, att, art, detDF, coeff_beta); + nodeApplyAscOrthoRadialGive(i_r, i_theta, grid_, DirBC_Interior_, smoother_color, x, rhs, temp, arr, att, art, + detDF, coeff_beta); } } diff --git a/tests/GMGPolar/solve_tests.cpp b/tests/GMGPolar/solve_tests.cpp index 7b9a9699..f7be2f40 100644 --- a/tests/GMGPolar/solve_tests.cpp +++ b/tests/GMGPolar/solve_tests.cpp @@ -510,6 +510,40 @@ using gmgpolar_test_suite = testing::Types< std::integral_constant, // expected_l2_error std::integral_constant, // expected_inf_error std::integral_constant // expected_residual_reduction + >, + std::tuple< + CzarnyGeometry, + SonnendruckerGyroCoefficients, + PolarR6_Boundary_CzarnyGeometry, + PolarR6_Sonnendrucker_CzarnyGeometry, + PolarR6_CzarnyGeometry, + std::integral_constant, // R0 + std::integral_constant, // Rmax + std::integral_constant, // nrExp + std::integral_constant, // nthetaExp + std::integral_constant, // refinementRadius + std::integral_constant, // anisotropicFactor + std::integral_constant, // divideBy2 + std::integral_constant, // verbose + std::integral_constant, // maxOpenMPThreads + std::integral_constant, // DirBC_Interior + std::integral_constant, // StencilDistributionMethod + std::integral_constant, // cacheDensityProfileCoefficient + std::integral_constant, // cacheDomainGeometry + std::integral_constant, // extrapolation + std::integral_constant, // maxLevels + std::integral_constant, // multigridCycle + std::integral_constant, // FMG + std::integral_constant, // FMG_iterations + std::integral_constant, // FMG_cycle + std::integral_constant, // maxIterations + std::integral_constant, // residualNormType + std::integral_constant, // absoluteTolerance + std::integral_constant, // relativeTolerance + std::integral_constant, // expected_iterations + std::integral_constant, // expected_l2_error + std::integral_constant, // expected_inf_error + std::integral_constant // expected_residual_reduction > >; // clang-format on