diff --git a/include/DirectSolver/DirectSolver-COO-MUMPS-Give/directSolverGive.h b/include/DirectSolver/DirectSolver-COO-MUMPS-Give/directSolverGive.h index 168a6eb7..b85cc98b 100644 --- a/include/DirectSolver/DirectSolver-COO-MUMPS-Give/directSolverGive.h +++ b/include/DirectSolver/DirectSolver-COO-MUMPS-Give/directSolverGive.h @@ -86,6 +86,10 @@ class DirectSolver_COO_MUMPS_Give : public DirectSolver // Retrieves the stencil for the solver matrix at the given radial index. const Stencil& getStencil(int i_r) const; + + void nodeBuildSolverMatrixGive(int i_r, int i_theta, const PolarGrid& grid, bool DirBC_Interior, + SparseMatrixCOO& solver_matrix, double arr, double att, double art, + double detDF, double coeff_beta); }; #endif \ No newline at end of file diff --git a/include/DirectSolver/DirectSolver-CSR-LU-Give/directSolverGiveCustomLU.h b/include/DirectSolver/DirectSolver-CSR-LU-Give/directSolverGiveCustomLU.h index bf4100c8..862d31f2 100644 --- a/include/DirectSolver/DirectSolver-CSR-LU-Give/directSolverGiveCustomLU.h +++ b/include/DirectSolver/DirectSolver-CSR-LU-Give/directSolverGiveCustomLU.h @@ -58,4 +58,8 @@ class DirectSolver_CSR_LU_Give : public DirectSolver const Stencil& getStencil(int i_r) const; int getStencilSize(int global_index) const; + + void nodeBuildSolverMatrixGive(int i_r, int i_theta, const PolarGrid& grid, const bool DirBC_Interior, + SparseMatrixCSR& solver_matrix, double arr, double att, double art, + double detDF, double coeff_beta); }; diff --git a/src/DirectSolver/DirectSolver-COO-MUMPS-Give/applySymmetryShift.cpp b/src/DirectSolver/DirectSolver-COO-MUMPS-Give/applySymmetryShift.cpp index 812f7964..02f5507b 100644 --- a/src/DirectSolver/DirectSolver-COO-MUMPS-Give/applySymmetryShift.cpp +++ b/src/DirectSolver/DirectSolver-COO-MUMPS-Give/applySymmetryShift.cpp @@ -122,35 +122,25 @@ void DirectSolver_COO_MUMPS_Give::applySymmetryShiftOuterBoundary(Vector } } -// clang-format off void DirectSolver_COO_MUMPS_Give::applySymmetryShift(Vector x) const { assert(std::ssize(x) == grid_.numberOfNodes()); assert(grid_.nr() >= 4); - if (num_omp_threads_ == 1) { - /* Single-threaded execution */ - if (DirBC_Interior_) { - applySymmetryShiftInnerBoundary(x); - } - applySymmetryShiftOuterBoundary(x); - } - else { - #pragma omp parallel sections num_threads(num_omp_threads_) + #pragma omp parallel sections num_threads(num_omp_threads_) + { + #pragma omp section { - #pragma omp section - { - if (DirBC_Interior_) { - applySymmetryShiftInnerBoundary(x); - } + if (DirBC_Interior_) { + applySymmetryShiftInnerBoundary(x); } + } - #pragma omp section - { - applySymmetryShiftOuterBoundary(x); - } + #pragma omp section + { + applySymmetryShiftOuterBoundary(x); } } } -// clang-format on + #endif diff --git a/src/DirectSolver/DirectSolver-COO-MUMPS-Give/buildSolverMatrix.cpp b/src/DirectSolver/DirectSolver-COO-MUMPS-Give/buildSolverMatrix.cpp index 63b8ce3b..54f58b64 100644 --- a/src/DirectSolver/DirectSolver-COO-MUMPS-Give/buildSolverMatrix.cpp +++ b/src/DirectSolver/DirectSolver-COO-MUMPS-Give/buildSolverMatrix.cpp @@ -4,766 +4,773 @@ #ifdef GMGPOLAR_USE_MUMPS - #define UPDATE_MATRIX_ELEMENT(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) - - #define NODE_BUILD_SOLVER_MATRIX_GIVE(i_r, i_theta, r, theta, grid, DirBC_Interior, solver_matrix, arr, att, art, \ - detDF, coeff_beta) \ - do { \ - int ptr, offset; \ - int row, col; \ - double val; \ - /* -------------------- */ \ - /* Node in the interior */ \ - /* -------------------- */ \ - if (i_r > 1 && i_r < grid.nr() - 2) { \ - const int i_theta_M1 = grid.wrapThetaIndex(i_theta - 1); \ - const int i_theta_P1 = grid.wrapThetaIndex(i_theta + 1); \ - \ - const double h1 = grid.radialSpacing(i_r - 1); \ - const double h2 = grid.radialSpacing(i_r); \ - const double k1 = grid.angularSpacing(i_theta_M1); \ - const double k2 = grid.angularSpacing(i_theta); \ - const double coeff1 = 0.5 * (k1 + k2) / h1; \ - const double coeff2 = 0.5 * (k1 + k2) / h2; \ - const double coeff3 = 0.5 * (h1 + h2) / k1; \ - const double coeff4 = 0.5 * (h1 + h2) / k2; \ - \ - const int center_nz_index = getSolverMatrixIndex(i_r, i_theta); \ - const int left_nz_index = getSolverMatrixIndex(i_r - 1, i_theta); \ - const int right_nz_index = getSolverMatrixIndex(i_r + 1, i_theta); \ - const int bottom_nz_index = getSolverMatrixIndex(i_r, i_theta_M1); \ - const int top_nz_index = getSolverMatrixIndex(i_r, i_theta_P1); \ - \ - const int center_index = grid.index(i_r, i_theta); \ - const int left_index = grid.index(i_r - 1, i_theta); \ - const int right_index = grid.index(i_r + 1, i_theta); \ - const int bottom_index = grid.index(i_r, i_theta_M1); \ - const int top_index = grid.index(i_r, i_theta_P1); \ - \ - /* Fill matrix row of (i,j) */ \ - row = center_index; \ - ptr = center_nz_index; \ - \ - const Stencil& CenterStencil = getStencil(i_r); \ - \ - offset = CenterStencil[StencilPosition::Center]; \ - col = center_index; \ - val = 0.25 * (h1 + h2) * (k1 + k2) * coeff_beta * fabs(detDF); /* beta_{i,j} */ \ - UPDATE_MATRIX_ELEMENT(solver_matrix, ptr, offset, row, col, val); \ - \ - offset = CenterStencil[StencilPosition::Left]; \ - col = left_index; \ - val = -coeff1 * arr; /* Left */ \ - UPDATE_MATRIX_ELEMENT(solver_matrix, ptr, offset, row, col, val); \ - \ - offset = CenterStencil[StencilPosition::Right]; \ - col = right_index; \ - val = -coeff2 * arr; /* Right */ \ - UPDATE_MATRIX_ELEMENT(solver_matrix, ptr, offset, row, col, val); \ - \ - offset = CenterStencil[StencilPosition::Bottom]; \ - col = bottom_index; \ - val = -coeff3 * att; /* Bottom */ \ - UPDATE_MATRIX_ELEMENT(solver_matrix, ptr, offset, row, col, val); \ - \ - offset = CenterStencil[StencilPosition::Top]; \ - col = top_index; \ - val = -coeff4 * att; /* Top */ \ - UPDATE_MATRIX_ELEMENT(solver_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) */ \ - UPDATE_MATRIX_ELEMENT(solver_matrix, ptr, offset, row, col, val); \ - \ - /* Fill matrix row of (i-1,j) */ \ - row = left_index; \ - ptr = left_nz_index; \ - \ - const Stencil& LeftStencil = getStencil(i_r - 1); \ - \ - offset = LeftStencil[StencilPosition::Right]; \ - col = center_index; \ - val = -coeff1 * arr; /* Right */ \ - UPDATE_MATRIX_ELEMENT(solver_matrix, ptr, offset, row, col, val); \ - \ - offset = LeftStencil[StencilPosition::Center]; \ - col = left_index; \ - val = +coeff1 * arr; /* Center: (Right) */ \ - UPDATE_MATRIX_ELEMENT(solver_matrix, ptr, offset, row, col, val); \ - \ - offset = LeftStencil[StencilPosition::TopRight]; \ - col = top_index; \ - val = -0.25 * art; /* Top Right */ \ - UPDATE_MATRIX_ELEMENT(solver_matrix, ptr, offset, row, col, val); \ - \ - offset = LeftStencil[StencilPosition::BottomRight]; \ - col = bottom_index; \ - val = +0.25 * art; /* Bottom Right */ \ - UPDATE_MATRIX_ELEMENT(solver_matrix, ptr, offset, row, col, val); \ - \ - /* Fill matrix row of (i+1,j) */ \ - row = right_index; \ - ptr = right_nz_index; \ - \ - const Stencil& RightStencil = getStencil(i_r + 1); \ - \ - offset = RightStencil[StencilPosition::Left]; \ - col = center_index; \ - val = -coeff2 * arr; /* Left */ \ - UPDATE_MATRIX_ELEMENT(solver_matrix, ptr, offset, row, col, val); \ - \ - offset = RightStencil[StencilPosition::Center]; \ - col = right_index; \ - val = +coeff2 * arr; /* Center: (Left) */ \ - UPDATE_MATRIX_ELEMENT(solver_matrix, ptr, offset, row, col, val); \ - \ - offset = RightStencil[StencilPosition::TopLeft]; \ - col = top_index; \ - val = +0.25 * art; /* Top Left */ \ - UPDATE_MATRIX_ELEMENT(solver_matrix, ptr, offset, row, col, val); \ - \ - offset = RightStencil[StencilPosition::BottomLeft]; \ - col = bottom_index; \ - val = -0.25 * art; /* Bottom Left */ \ - UPDATE_MATRIX_ELEMENT(solver_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::Top]; \ - col = center_index; \ - val = -coeff3 * att; /* Top */ \ - UPDATE_MATRIX_ELEMENT(solver_matrix, ptr, offset, row, col, val); \ - \ - offset = BottomStencil[StencilPosition::Center]; \ - col = bottom_index; \ - val = +coeff3 * att; /* Center: (Top) */ \ - UPDATE_MATRIX_ELEMENT(solver_matrix, ptr, offset, row, col, val); \ - \ - offset = BottomStencil[StencilPosition::TopRight]; \ - col = right_index; \ - val = -0.25 * art; /* Top Right */ \ - UPDATE_MATRIX_ELEMENT(solver_matrix, ptr, offset, row, col, val); \ - \ - offset = BottomStencil[StencilPosition::TopLeft]; \ - col = left_index; \ - val = +0.25 * art; /* Top Left */ \ - UPDATE_MATRIX_ELEMENT(solver_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::Bottom]; \ - col = center_index; \ - val = -coeff4 * att; /* Bottom */ \ - UPDATE_MATRIX_ELEMENT(solver_matrix, ptr, offset, row, col, val); \ - \ - offset = TopStencil[StencilPosition::Center]; \ - col = top_index; \ - val = +coeff4 * att; /* Center: (Bottom) */ \ - UPDATE_MATRIX_ELEMENT(solver_matrix, ptr, offset, row, col, val); \ - \ - offset = TopStencil[StencilPosition::BottomRight]; \ - col = right_index; \ - val = +0.25 * art; /* Bottom Right */ \ - UPDATE_MATRIX_ELEMENT(solver_matrix, ptr, offset, row, col, val); \ - \ - offset = TopStencil[StencilPosition::BottomLeft]; \ - col = left_index; \ - val = -0.25 * art; /* Bottom Left */ \ - UPDATE_MATRIX_ELEMENT(solver_matrix, ptr, offset, row, col, val); \ - } \ - /* -------------------------- */ \ - /* Node on the inner boundary */ \ - /* -------------------------- */ \ - else if (i_r == 0) { \ - /* ------------------------------------------------ */ \ - /* Case 1: Dirichlet boundary on the inner boundary */ \ - /* ------------------------------------------------ */ \ - if (DirBC_Interior) { \ - const double h2 = grid.radialSpacing(i_r); \ - const double k1 = grid.angularSpacing(i_theta - 1); \ - const double k2 = grid.angularSpacing(i_theta); \ - const double coeff2 = 0.5 * (k1 + k2) / h2; \ - \ - const int i_theta_M1 = grid.wrapThetaIndex(i_theta - 1); \ - const int i_theta_P1 = grid.wrapThetaIndex(i_theta + 1); \ - \ - const int center_nz_index = getSolverMatrixIndex(i_r, i_theta); \ - const int right_nz_index = getSolverMatrixIndex(i_r + 1, i_theta); \ - \ - const int center_index = grid.index(i_r, i_theta); \ - const int right_index = grid.index(i_r + 1, i_theta); \ - const int bottom_index = grid.index(i_r, i_theta_M1); \ - const int top_index = grid.index(i_r, i_theta_P1); \ - \ - /* Fill matrix row of (i,j) */ \ - row = center_index; \ - ptr = center_nz_index; \ - \ - const Stencil& CenterStencil = getStencil(i_r); \ - \ - offset = CenterStencil[StencilPosition::Center]; \ - col = center_index; \ - val = 1.0; \ - UPDATE_MATRIX_ELEMENT(solver_matrix, ptr, offset, row, col, val); \ - \ - /* Fill matrix row of (i+1,j) */ \ - row = right_index; \ - ptr = right_nz_index; \ - \ - const Stencil& RightStencil = getStencil(i_r + 1); \ - \ - /* Left REMOVED: Moved to the right hand side to make the matrix symmetric */ \ - \ - offset = RightStencil[StencilPosition::Center]; \ - col = right_index; \ - val = +coeff2 * arr; /* Center: (Left) */ \ - UPDATE_MATRIX_ELEMENT(solver_matrix, ptr, offset, row, col, val); \ - \ - /* TopLeft REMOVED: Moved to the right hand side to make the matrix symmetric */ \ - \ - /* BottomLeft REMOVED: Moved to the right hand side to make the matrix symmetric */ \ - } \ - 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. */ \ - const int i_theta_M1 = grid.wrapThetaIndex(i_theta - 1); \ - const int i_theta_P1 = grid.wrapThetaIndex(i_theta + 1); \ - \ - assert(grid_.ntheta() % 2 == 0); \ - const int i_theta_AcrossOrigin = grid.wrapThetaIndex(i_theta + grid.ntheta() / 2); \ - \ - double h1 = 2.0 * grid.radius(0); \ - double h2 = grid.radialSpacing(i_r); \ - double k1 = grid.angularSpacing(i_theta_M1); \ - 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 center_nz_index = getSolverMatrixIndex(i_r, i_theta); \ - const int left_nz_index = getSolverMatrixIndex(i_r, i_theta_AcrossOrigin); \ - const int right_nz_index = getSolverMatrixIndex(i_r + 1, i_theta); \ - const int bottom_nz_index = getSolverMatrixIndex(i_r, i_theta_M1); \ - const int top_nz_index = getSolverMatrixIndex(i_r, i_theta_P1); \ - \ - const int center_index = grid.index(i_r, i_theta); \ - const int left_index = grid.index(i_r, i_theta_AcrossOrigin); \ - const int right_index = grid.index(i_r + 1, i_theta); \ - const int bottom_index = grid.index(i_r, i_theta_M1); \ - const int top_index = grid.index(i_r, i_theta_P1); \ - \ - /* Fill matrix row of (i,j) */ \ - row = center_index; \ - ptr = center_nz_index; \ - const Stencil& CenterStencil = getStencil(i_r); \ - \ - offset = CenterStencil[StencilPosition::Center]; \ - col = center_index; \ - val = 0.25 * (h1 + h2) * (k1 + k2) * coeff_beta * fabs(detDF); /* beta_{i,j} */ \ - UPDATE_MATRIX_ELEMENT(solver_matrix, ptr, offset, row, col, val); \ - \ - offset = CenterStencil[StencilPosition::Left]; \ - col = left_index; \ - val = -coeff1 * arr; /* Left */ \ - UPDATE_MATRIX_ELEMENT(solver_matrix, ptr, offset, row, col, val); \ - \ - offset = CenterStencil[StencilPosition::Right]; \ - col = right_index; \ - val = -coeff2 * arr; /* Right */ \ - UPDATE_MATRIX_ELEMENT(solver_matrix, ptr, offset, row, col, val); \ - \ - offset = CenterStencil[StencilPosition::Bottom]; \ - col = bottom_index; \ - val = -coeff3 * att; /* Bottom */ \ - UPDATE_MATRIX_ELEMENT(solver_matrix, ptr, offset, row, col, val); \ - \ - offset = CenterStencil[StencilPosition::Top]; \ - col = top_index; \ - val = -coeff4 * att; /* Top */ \ - UPDATE_MATRIX_ELEMENT(solver_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) */ \ - UPDATE_MATRIX_ELEMENT(solver_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*/ \ - UPDATE_MATRIX_ELEMENT(solver_matrix, ptr, offset, row, col, val); \ - \ - offset = LeftStencil[StencilPosition::Center]; \ - col = left_index; \ - val = +coeff1 * arr; /* Center: (Right) -> Center: (Left) */ \ - UPDATE_MATRIX_ELEMENT(solver_matrix, ptr, offset, row, col, val); \ - \ - /* Top Right -> Bottom Left: REMOVED DUE TO ARTIFICAL 7 POINT STENCIL */ \ - \ - /* Bottom Right -> Top Left: REMOVED DUE TO ARTIFICAL 7 POINT STENCIL */ \ - \ - /* Fill matrix row of (i+1,j) */ \ - row = right_index; \ - ptr = right_nz_index; \ - const Stencil& RightStencil = getStencil(i_r + 1); \ - \ - offset = RightStencil[StencilPosition::Left]; \ - col = center_index; \ - val = -coeff2 * arr; /* Left */ \ - UPDATE_MATRIX_ELEMENT(solver_matrix, ptr, offset, row, col, val); \ - \ - offset = RightStencil[StencilPosition::Center]; \ - col = right_index; \ - val = +coeff2 * arr; /* Center: (Left) */ \ - UPDATE_MATRIX_ELEMENT(solver_matrix, ptr, offset, row, col, val); \ - \ - offset = RightStencil[StencilPosition::TopLeft]; \ - col = top_index; \ - val = +0.25 * art; /* Top Left */ \ - UPDATE_MATRIX_ELEMENT(solver_matrix, ptr, offset, row, col, val); \ - \ - offset = RightStencil[StencilPosition::BottomLeft]; \ - col = bottom_index; \ - val = -0.25 * art; /* Bottom Left */ \ - UPDATE_MATRIX_ELEMENT(solver_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::Top]; \ - col = center_index; \ - val = -coeff3 * att; /* Top */ \ - UPDATE_MATRIX_ELEMENT(solver_matrix, ptr, offset, row, col, val); \ - \ - offset = BottomStencil[StencilPosition::Center]; \ - col = bottom_index; \ - val = +coeff3 * att; /* Center: (Top) */ \ - UPDATE_MATRIX_ELEMENT(solver_matrix, ptr, offset, row, col, val); \ - \ - offset = BottomStencil[StencilPosition::TopRight]; \ - col = right_index; \ - val = -0.25 * art; /* Top Right */ \ - UPDATE_MATRIX_ELEMENT(solver_matrix, ptr, offset, row, col, val); \ - \ - /* TopLeft REMOVED DUE TO ARTIFICAL 7 POINT STENCIL */ \ - \ - /* Fill matrix row of (i,j+1) */ \ - row = top_index; \ - ptr = top_nz_index; \ - \ - const Stencil& TopStencil = CenterStencil; \ - \ - offset = TopStencil[StencilPosition::Bottom]; \ - col = center_index; \ - val = -coeff4 * att; /* Bottom */ \ - UPDATE_MATRIX_ELEMENT(solver_matrix, ptr, offset, row, col, val); \ - \ - offset = TopStencil[StencilPosition::Center]; \ - col = top_index; \ - val = +coeff4 * att; /* Center: (Bottom) */ \ - UPDATE_MATRIX_ELEMENT(solver_matrix, ptr, offset, row, col, val); \ - \ - offset = TopStencil[StencilPosition::BottomRight]; \ - col = right_index; \ - val = +0.25 * art; /* Bottom Right */ \ - UPDATE_MATRIX_ELEMENT(solver_matrix, ptr, offset, row, col, val); \ - \ - /* BottomLeft REMOVED DUE TO ARTIFICAL 7 POINT STENCIL */ \ - } \ - } \ - /* ------------------------------- */ \ - /* Node next to the inner boundary */ \ - /* ------------------------------- */ \ - else if (i_r == 1) { \ - const double h1 = grid.radialSpacing(i_r - 1); \ - const double h2 = grid.radialSpacing(i_r); \ - const double k1 = grid.angularSpacing(i_theta - 1); \ - const double k2 = grid.angularSpacing(i_theta); \ - const double coeff1 = 0.5 * (k1 + k2) / h1; \ - const double coeff2 = 0.5 * (k1 + k2) / h2; \ - const double coeff3 = 0.5 * (h1 + h2) / k1; \ - const 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_nz_index = getSolverMatrixIndex(i_r, i_theta); \ - const int left_nz_index = getSolverMatrixIndex(i_r - 1, i_theta); \ - const int right_nz_index = getSolverMatrixIndex(i_r + 1, i_theta); \ - const int bottom_nz_index = getSolverMatrixIndex(i_r, i_theta_M1); \ - const int top_nz_index = getSolverMatrixIndex(i_r, i_theta_P1); \ - \ - const int center_index = grid.index(i_r, i_theta); \ - const int left_index = grid.index(i_r - 1, i_theta); \ - const int right_index = grid.index(i_r + 1, i_theta); \ - const int bottom_index = grid.index(i_r, i_theta_M1); \ - const int top_index = grid.index(i_r, i_theta_P1); \ - \ - /* Fill matrix row of (i,j) */ \ - row = center_index; \ - ptr = center_nz_index; \ - \ - const Stencil& CenterStencil = getStencil(i_r); \ - \ - offset = CenterStencil[StencilPosition::Center]; \ - col = center_index; \ - val = 0.25 * (h1 + h2) * (k1 + k2) * coeff_beta * fabs(detDF); /* beta_{i,j} */ \ - UPDATE_MATRIX_ELEMENT(solver_matrix, ptr, offset, row, col, val); \ - \ - /* REMOVED: Moved to the right hand side to make the matrix symmetric */ \ - if (!DirBC_Interior) { \ - offset = CenterStencil[StencilPosition::Left]; \ - col = left_index; \ - val = -coeff1 * arr; /* Left */ \ - UPDATE_MATRIX_ELEMENT(solver_matrix, ptr, offset, row, col, val); \ - } \ - \ - offset = CenterStencil[StencilPosition::Right]; \ - col = right_index; \ - val = -coeff2 * arr; /* Right */ \ - UPDATE_MATRIX_ELEMENT(solver_matrix, ptr, offset, row, col, val); \ - \ - offset = CenterStencil[StencilPosition::Bottom]; \ - col = bottom_index; \ - val = -coeff3 * att; /* Bottom */ \ - UPDATE_MATRIX_ELEMENT(solver_matrix, ptr, offset, row, col, val); \ - \ - offset = CenterStencil[StencilPosition::Top]; \ - col = top_index; \ - val = -coeff4 * att; /* Top */ \ - UPDATE_MATRIX_ELEMENT(solver_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) */ \ - UPDATE_MATRIX_ELEMENT(solver_matrix, ptr, offset, row, col, val); \ - \ - if (!DirBC_Interior) { /* Don't give to the inner Dirichlet boundary! */ \ - /* Fill matrix row of (i-1,j) */ \ - row = left_index; \ - ptr = left_nz_index; \ - \ - const Stencil& LeftStencil = getStencil(i_r - 1); \ - \ - offset = LeftStencil[StencilPosition::Right]; \ - col = center_index; \ - val = -coeff1 * arr; /* Right */ \ - UPDATE_MATRIX_ELEMENT(solver_matrix, ptr, offset, row, col, val); \ - \ - offset = LeftStencil[StencilPosition::Center]; \ - col = left_index; \ - val = +coeff1 * arr; /* Center: (Right) */ \ - UPDATE_MATRIX_ELEMENT(solver_matrix, ptr, offset, row, col, val); \ - \ - offset = LeftStencil[StencilPosition::TopRight]; \ - col = top_index; \ - val = -0.25 * art; /* Top Right */ \ - UPDATE_MATRIX_ELEMENT(solver_matrix, ptr, offset, row, col, val); \ - \ - offset = LeftStencil[StencilPosition::BottomRight]; \ - col = bottom_index; \ - val = +0.25 * art; /* Bottom Right */ \ - UPDATE_MATRIX_ELEMENT(solver_matrix, ptr, offset, row, col, val); \ - } \ - \ - /* Fill matrix row of (i+1,j) */ \ - row = right_index; \ - ptr = right_nz_index; \ - \ - const Stencil& RightStencil = getStencil(i_r + 1); \ - \ - offset = RightStencil[StencilPosition::Left]; \ - col = center_index; \ - val = -coeff2 * arr; /* Left */ \ - UPDATE_MATRIX_ELEMENT(solver_matrix, ptr, offset, row, col, val); \ - \ - offset = RightStencil[StencilPosition::Center]; \ - col = right_index; \ - val = +coeff2 * arr; /* Center: (Left) */ \ - UPDATE_MATRIX_ELEMENT(solver_matrix, ptr, offset, row, col, val); \ - \ - offset = RightStencil[StencilPosition::TopLeft]; \ - col = top_index; \ - val = +0.25 * art; /* Top Left */ \ - UPDATE_MATRIX_ELEMENT(solver_matrix, ptr, offset, row, col, val); \ - \ - offset = RightStencil[StencilPosition::BottomLeft]; \ - col = bottom_index; \ - val = -0.25 * art; /* Bottom Left */ \ - UPDATE_MATRIX_ELEMENT(solver_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::Top]; \ - col = center_index; \ - val = -coeff3 * att; /* Top */ \ - UPDATE_MATRIX_ELEMENT(solver_matrix, ptr, offset, row, col, val); \ - \ - offset = BottomStencil[StencilPosition::Center]; \ - col = bottom_index; \ - val = +coeff3 * att; /* Center: (Top) */ \ - UPDATE_MATRIX_ELEMENT(solver_matrix, ptr, offset, row, col, val); \ - \ - offset = BottomStencil[StencilPosition::TopRight]; \ - col = right_index; \ - val = -0.25 * art; /* Top Right */ \ - UPDATE_MATRIX_ELEMENT(solver_matrix, ptr, offset, row, col, val); \ - \ - /* REMOVED: Moved to the right hand side to make the matrix symmetric */ \ - if (!DirBC_Interior) { \ - offset = BottomStencil[StencilPosition::TopLeft]; \ - col = left_index; \ - val = +0.25 * art; /* Top Left */ \ - UPDATE_MATRIX_ELEMENT(solver_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::Bottom]; \ - col = center_index; \ - val = -coeff4 * att; /* Bottom */ \ - UPDATE_MATRIX_ELEMENT(solver_matrix, ptr, offset, row, col, val); \ - \ - offset = TopStencil[StencilPosition::Center]; \ - col = top_index; \ - val = +coeff4 * att; /* Center: (Bottom) */ \ - UPDATE_MATRIX_ELEMENT(solver_matrix, ptr, offset, row, col, val); \ - \ - offset = TopStencil[StencilPosition::BottomRight]; \ - col = right_index; \ - val = +0.25 * art; /* Bottom Right */ \ - UPDATE_MATRIX_ELEMENT(solver_matrix, ptr, offset, row, col, val); \ - \ - /* REMOVED: Moved to the right hand side to make the matrix symmetric */ \ - if (!DirBC_Interior) { \ - offset = TopStencil[StencilPosition::BottomLeft]; \ - col = left_index; \ - val = -0.25 * art; /* Bottom Left */ \ - UPDATE_MATRIX_ELEMENT(solver_matrix, ptr, offset, row, col, val); \ - } \ - } \ - /* ------------------------------- */ \ - /* Node next to the outer boundary */ \ - /* ------------------------------- */ \ - else if (i_r == grid.nr() - 2) { \ - const double h1 = grid.radialSpacing(i_r - 1); \ - const double h2 = grid.radialSpacing(i_r); \ - const double k1 = grid.angularSpacing(i_theta - 1); \ - const double k2 = grid.angularSpacing(i_theta); \ - const double coeff1 = 0.5 * (k1 + k2) / h1; \ - const double coeff2 = 0.5 * (k1 + k2) / h2; \ - const double coeff3 = 0.5 * (h1 + h2) / k1; \ - const 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_nz_index = getSolverMatrixIndex(i_r, i_theta); \ - const int left_nz_index = getSolverMatrixIndex(i_r - 1, i_theta); \ - const int right_nz_index = getSolverMatrixIndex(i_r + 1, i_theta); \ - const int bottom_nz_index = getSolverMatrixIndex(i_r, i_theta_M1); \ - const int top_nz_index = getSolverMatrixIndex(i_r, i_theta_P1); \ - \ - const int center_index = grid.index(i_r, i_theta); \ - const int left_index = grid.index(i_r - 1, i_theta); \ - const int right_index = grid.index(i_r + 1, i_theta); \ - const int bottom_index = grid.index(i_r, i_theta_M1); \ - const int top_index = grid.index(i_r, i_theta_P1); \ - \ - /* Fill matrix row of (i,j) */ \ - row = center_index; \ - ptr = center_nz_index; \ - \ - const Stencil& CenterStencil = getStencil(i_r); \ - \ - offset = CenterStencil[StencilPosition::Center]; \ - col = center_index; \ - val = 0.25 * (h1 + h2) * (k1 + k2) * coeff_beta * fabs(detDF); /* beta_{i,j} */ \ - UPDATE_MATRIX_ELEMENT(solver_matrix, ptr, offset, row, col, val); \ - \ - offset = CenterStencil[StencilPosition::Left]; \ - col = left_index; \ - val = -coeff1 * arr; /* Left */ \ - UPDATE_MATRIX_ELEMENT(solver_matrix, ptr, offset, row, col, val); \ - \ - /* Right REMOVED: Moved to the right hand side to make the matrix symmetric */ \ - \ - offset = CenterStencil[StencilPosition::Bottom]; \ - col = bottom_index; \ - val = -coeff3 * att; /* Bottom */ \ - UPDATE_MATRIX_ELEMENT(solver_matrix, ptr, offset, row, col, val); \ - \ - offset = CenterStencil[StencilPosition::Top]; \ - col = top_index; \ - val = -coeff4 * att; /* Top */ \ - UPDATE_MATRIX_ELEMENT(solver_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) */ \ - UPDATE_MATRIX_ELEMENT(solver_matrix, ptr, offset, row, col, val); \ - \ - /* Fill matrix row of (i-1,j) */ \ - row = left_index; \ - ptr = left_nz_index; \ - \ - const Stencil& LeftStencil = getStencil(i_r - 1); \ - \ - offset = LeftStencil[StencilPosition::Right]; \ - col = center_index; \ - val = -coeff1 * arr; /* Right */ \ - UPDATE_MATRIX_ELEMENT(solver_matrix, ptr, offset, row, col, val); \ - \ - offset = LeftStencil[StencilPosition::Center]; \ - col = left_index; \ - val = coeff1 * arr; /* Center: (Right) */ \ - UPDATE_MATRIX_ELEMENT(solver_matrix, ptr, offset, row, col, val); \ - \ - offset = LeftStencil[StencilPosition::TopRight]; \ - col = top_index; \ - val = -0.25 * art; /* Top Right */ \ - UPDATE_MATRIX_ELEMENT(solver_matrix, ptr, offset, row, col, val); \ - \ - offset = LeftStencil[StencilPosition::BottomRight]; \ - col = bottom_index; \ - val = 0.25 * art; /* Bottom Right */ \ - UPDATE_MATRIX_ELEMENT(solver_matrix, ptr, offset, row, col, val); \ - \ - /* Fill matrix row of (i+1,j) */ \ - /* Don't give to the outer dirichlet boundary! */ \ - \ - /* Fill matrix row of (i,j-1) */ \ - row = bottom_index; \ - ptr = bottom_nz_index; \ - \ - const Stencil& BottomStencil = CenterStencil; \ - \ - offset = BottomStencil[StencilPosition::Top]; \ - col = center_index; \ - val = -coeff3 * att; /* Top */ \ - UPDATE_MATRIX_ELEMENT(solver_matrix, ptr, offset, row, col, val); \ - \ - offset = BottomStencil[StencilPosition::Center]; \ - col = bottom_index; \ - val = coeff3 * att; /* Center: (Top) */ \ - UPDATE_MATRIX_ELEMENT(solver_matrix, ptr, offset, row, col, val); \ - \ - /* TopRight REMOVED: Moved to the right hand side to make the matrix symmetric */ \ - \ - offset = BottomStencil[StencilPosition::TopLeft]; \ - col = left_index; \ - val = 0.25 * art; /* Top Left */ \ - UPDATE_MATRIX_ELEMENT(solver_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::Bottom]; \ - col = center_index; \ - val = -coeff4 * att; /* Bottom */ \ - UPDATE_MATRIX_ELEMENT(solver_matrix, ptr, offset, row, col, val); \ - \ - offset = TopStencil[StencilPosition::Center]; \ - col = top_index; \ - val = coeff4 * att; /* Center: (Bottom) */ \ - UPDATE_MATRIX_ELEMENT(solver_matrix, ptr, offset, row, col, val); \ - \ - /* BottomRight REMOVED: Moved to the right hand side to make the matrix symmetric */ \ - \ - offset = TopStencil[StencilPosition::BottomLeft]; \ - col = left_index; \ - val = -0.25 * art; /* Bottom Left */ \ - UPDATE_MATRIX_ELEMENT(solver_matrix, ptr, offset, row, col, val); \ - } \ - /* ------------------------------------ */ \ - /* Node on the outer dirichlet boundary */ \ - /* ------------------------------------ */ \ - else if (i_r == grid.nr() - 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; \ - \ - const int i_theta_M1 = grid.wrapThetaIndex(i_theta - 1); \ - const int i_theta_P1 = grid.wrapThetaIndex(i_theta + 1); \ - \ - const int center_nz_index = getSolverMatrixIndex(i_r, i_theta); \ - const int left_nz_index = getSolverMatrixIndex(i_r - 1, i_theta); \ - \ - const int center_index = grid.index(i_r, i_theta); \ - const int left_index = grid.index(i_r - 1, i_theta); \ - const int bottom_index = grid.index(i_r, i_theta_M1); \ - const int top_index = grid.index(i_r, i_theta_P1); \ - \ - /* Fill matrix row of (i,j) */ \ - row = center_index; \ - ptr = center_nz_index; \ - \ - const Stencil& CenterStencil = getStencil(i_r); \ - \ - offset = CenterStencil[StencilPosition::Center]; \ - col = center_index; \ - val = 1.0; \ - UPDATE_MATRIX_ELEMENT(solver_matrix, ptr, offset, row, col, val); \ - \ - /* Give value to the interior nodes! */ \ - /* Fill matrix row of (i-1,j) */ \ - row = left_index; \ - ptr = left_nz_index; \ - const Stencil& LeftStencil = getStencil(i_r - 1); \ - \ - /* Right REMOVED: Moved to the right hand side to make the matrix symmetric */ \ - \ - offset = LeftStencil[StencilPosition::Center]; \ - col = left_index; \ - val = coeff1 * arr; /* Center: (Right) */ \ - UPDATE_MATRIX_ELEMENT(solver_matrix, ptr, offset, row, col, val); \ - \ - /* TopRight REMOVED: Moved to the right hand side to make the matrix symmetric */ \ - \ - /* BottomRight REMOVED: Moved to the right hand side to make the matrix symmetric */ \ - } \ - } while (0) +inline void updateMatrixElement(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; +} + +void DirectSolver_COO_MUMPS_Give::nodeBuildSolverMatrixGive(int i_r, int i_theta, const PolarGrid& grid, + bool DirBC_Interior, SparseMatrixCOO& solver_matrix, + double arr, double att, double art, double detDF, + double coeff_beta) +{ + int ptr, offset; + int row, col; + double val; + /* -------------------- */ + /* Node in the interior */ + /* -------------------- */ + if (i_r > 1 && i_r < grid.nr() - 2) { + const int i_theta_M1 = grid.wrapThetaIndex(i_theta - 1); + const int i_theta_P1 = grid.wrapThetaIndex(i_theta + 1); + + const double h1 = grid.radialSpacing(i_r - 1); + const double h2 = grid.radialSpacing(i_r); + const double k1 = grid.angularSpacing(i_theta_M1); + const double k2 = grid.angularSpacing(i_theta); + + const double coeff1 = 0.5 * (k1 + k2) / h1; + const double coeff2 = 0.5 * (k1 + k2) / h2; + const double coeff3 = 0.5 * (h1 + h2) / k1; + const double coeff4 = 0.5 * (h1 + h2) / k2; + + const int center_nz_index = getSolverMatrixIndex(i_r, i_theta); + const int left_nz_index = getSolverMatrixIndex(i_r - 1, i_theta); + const int right_nz_index = getSolverMatrixIndex(i_r + 1, i_theta); + const int bottom_nz_index = getSolverMatrixIndex(i_r, i_theta_M1); + const int top_nz_index = getSolverMatrixIndex(i_r, i_theta_P1); + + const int center_index = grid.index(i_r, i_theta); + const int left_index = grid.index(i_r - 1, i_theta); + const int right_index = grid.index(i_r + 1, i_theta); + const int bottom_index = grid.index(i_r, i_theta_M1); + const int top_index = grid.index(i_r, i_theta_P1); + + /* Fill matrix row of (i,j) */ + row = center_index; + ptr = center_nz_index; + + const Stencil& CenterStencil = getStencil(i_r); + + offset = CenterStencil[StencilPosition::Center]; + col = center_index; + val = 0.25 * (h1 + h2) * (k1 + k2) * coeff_beta * std::fabs(detDF); /* beta_{i,j} */ + updateMatrixElement(solver_matrix, ptr, offset, row, col, val); + + offset = CenterStencil[StencilPosition::Left]; + col = left_index; + val = -coeff1 * arr; /* Left */ + updateMatrixElement(solver_matrix, ptr, offset, row, col, val); + + offset = CenterStencil[StencilPosition::Right]; + col = right_index; + val = -coeff2 * arr; /* Right */ + updateMatrixElement(solver_matrix, ptr, offset, row, col, val); + + offset = CenterStencil[StencilPosition::Bottom]; + col = bottom_index; + val = -coeff3 * att; /* Bottom */ + updateMatrixElement(solver_matrix, ptr, offset, row, col, val); + + offset = CenterStencil[StencilPosition::Top]; + col = top_index; + val = -coeff4 * att; /* Top */ + updateMatrixElement(solver_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) */ + updateMatrixElement(solver_matrix, ptr, offset, row, col, val); + + /* Fill matrix row of (i-1,j) */ + row = left_index; + ptr = left_nz_index; + + const Stencil& LeftStencil = getStencil(i_r - 1); + + offset = LeftStencil[StencilPosition::Right]; + col = center_index; + val = -coeff1 * arr; /* Right */ + updateMatrixElement(solver_matrix, ptr, offset, row, col, val); + + offset = LeftStencil[StencilPosition::Center]; + col = left_index; + val = +coeff1 * arr; /* Center: (Right) */ + updateMatrixElement(solver_matrix, ptr, offset, row, col, val); + + offset = LeftStencil[StencilPosition::TopRight]; + col = top_index; + val = -0.25 * art; /* Top Right */ + updateMatrixElement(solver_matrix, ptr, offset, row, col, val); + + offset = LeftStencil[StencilPosition::BottomRight]; + col = bottom_index; + val = +0.25 * art; /* Bottom Right */ + updateMatrixElement(solver_matrix, ptr, offset, row, col, val); + + /* Fill matrix row of (i+1,j) */ + row = right_index; + ptr = right_nz_index; + + const Stencil& RightStencil = getStencil(i_r + 1); + + offset = RightStencil[StencilPosition::Left]; + col = center_index; + val = -coeff2 * arr; /* Left */ + updateMatrixElement(solver_matrix, ptr, offset, row, col, val); + + offset = RightStencil[StencilPosition::Center]; + col = right_index; + val = +coeff2 * arr; /* Center: (Left) */ + updateMatrixElement(solver_matrix, ptr, offset, row, col, val); + + offset = RightStencil[StencilPosition::TopLeft]; + col = top_index; + val = +0.25 * art; /* Top Left */ + updateMatrixElement(solver_matrix, ptr, offset, row, col, val); + + offset = RightStencil[StencilPosition::BottomLeft]; + col = bottom_index; + val = -0.25 * art; /* Bottom Left */ + updateMatrixElement(solver_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::Top]; + col = center_index; + val = -coeff3 * att; /* Top */ + updateMatrixElement(solver_matrix, ptr, offset, row, col, val); + + offset = BottomStencil[StencilPosition::Center]; + col = bottom_index; + val = +coeff3 * att; /* Center: (Top) */ + updateMatrixElement(solver_matrix, ptr, offset, row, col, val); + + offset = BottomStencil[StencilPosition::TopRight]; + col = right_index; + val = -0.25 * art; /* Top Right */ + updateMatrixElement(solver_matrix, ptr, offset, row, col, val); + + offset = BottomStencil[StencilPosition::TopLeft]; + col = left_index; + val = +0.25 * art; /* Top Left */ + updateMatrixElement(solver_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::Bottom]; + col = center_index; + val = -coeff4 * att; /* Bottom */ + updateMatrixElement(solver_matrix, ptr, offset, row, col, val); + + offset = TopStencil[StencilPosition::Center]; + col = top_index; + val = +coeff4 * att; /* Center: (Bottom) */ + updateMatrixElement(solver_matrix, ptr, offset, row, col, val); + + offset = TopStencil[StencilPosition::BottomRight]; + col = right_index; + val = +0.25 * art; /* Bottom Right */ + updateMatrixElement(solver_matrix, ptr, offset, row, col, val); + + offset = TopStencil[StencilPosition::BottomLeft]; + col = left_index; + val = -0.25 * art; /* Bottom Left */ + updateMatrixElement(solver_matrix, ptr, offset, row, col, val); + } + /* -------------------------- */ + /* Node on the inner boundary */ + /* -------------------------- */ + else if (i_r == 0) { + /* ------------------------------------------------ */ + /* Case 1: Dirichlet boundary on the inner boundary */ + /* ------------------------------------------------ */ + if (DirBC_Interior) { + const double h2 = grid.radialSpacing(i_r); + const double k1 = grid.angularSpacing(i_theta - 1); + const double k2 = grid.angularSpacing(i_theta); + const double coeff2 = 0.5 * (k1 + k2) / h2; + + const int i_theta_M1 = grid.wrapThetaIndex(i_theta - 1); + const int i_theta_P1 = grid.wrapThetaIndex(i_theta + 1); + + const int center_nz_index = getSolverMatrixIndex(i_r, i_theta); + const int right_nz_index = getSolverMatrixIndex(i_r + 1, i_theta); + + const int center_index = grid.index(i_r, i_theta); + const int right_index = grid.index(i_r + 1, i_theta); + const int bottom_index = grid.index(i_r, i_theta_M1); + const int top_index = grid.index(i_r, i_theta_P1); + + /* Fill matrix row of (i,j) */ + row = center_index; + ptr = center_nz_index; + + const Stencil& CenterStencil = getStencil(i_r); + + offset = CenterStencil[StencilPosition::Center]; + col = center_index; + val = 1.0; + updateMatrixElement(solver_matrix, ptr, offset, row, col, val); + + /* Fill matrix row of (i+1,j) */ + row = right_index; + ptr = right_nz_index; + + const Stencil& RightStencil = getStencil(i_r + 1); + + /* Left REMOVED: Moved to the right hand side to make the matrix symmetric */ + + offset = RightStencil[StencilPosition::Center]; + col = right_index; + val = +coeff2 * arr; /* Center: (Left) */ + updateMatrixElement(solver_matrix, ptr, offset, row, col, val); + + /* TopLeft REMOVED: Moved to the right hand side to make the matrix symmetric */ + + /* BottomLeft REMOVED: Moved to the right hand side to make the matrix symmetric */ + } + 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. */ + const int i_theta_M1 = grid.wrapThetaIndex(i_theta - 1); + const int i_theta_P1 = grid.wrapThetaIndex(i_theta + 1); + + assert(grid_.ntheta() % 2 == 0); + const int i_theta_AcrossOrigin = grid.wrapThetaIndex(i_theta + grid.ntheta() / 2); + + double h1 = 2.0 * grid.radius(0); + double h2 = grid.radialSpacing(i_r); + double k1 = grid.angularSpacing(i_theta_M1); + 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 center_nz_index = getSolverMatrixIndex(i_r, i_theta); + const int left_nz_index = getSolverMatrixIndex(i_r, i_theta_AcrossOrigin); + const int right_nz_index = getSolverMatrixIndex(i_r + 1, i_theta); + const int bottom_nz_index = getSolverMatrixIndex(i_r, i_theta_M1); + const int top_nz_index = getSolverMatrixIndex(i_r, i_theta_P1); + + const int center_index = grid.index(i_r, i_theta); + const int left_index = grid.index(i_r, i_theta_AcrossOrigin); + const int right_index = grid.index(i_r + 1, i_theta); + const int bottom_index = grid.index(i_r, i_theta_M1); + const int top_index = grid.index(i_r, i_theta_P1); + + /* Fill matrix row of (i,j) */ + row = center_index; + ptr = center_nz_index; + const Stencil& CenterStencil = getStencil(i_r); + + offset = CenterStencil[StencilPosition::Center]; + col = center_index; + val = 0.25 * (h1 + h2) * (k1 + k2) * coeff_beta * std::fabs(detDF); /* beta_{i,j} */ + updateMatrixElement(solver_matrix, ptr, offset, row, col, val); + + offset = CenterStencil[StencilPosition::Left]; + col = left_index; + val = -coeff1 * arr; /* Left */ + updateMatrixElement(solver_matrix, ptr, offset, row, col, val); + + offset = CenterStencil[StencilPosition::Right]; + col = right_index; + val = -coeff2 * arr; /* Right */ + updateMatrixElement(solver_matrix, ptr, offset, row, col, val); + + offset = CenterStencil[StencilPosition::Bottom]; + col = bottom_index; + val = -coeff3 * att; /* Bottom */ + updateMatrixElement(solver_matrix, ptr, offset, row, col, val); + + offset = CenterStencil[StencilPosition::Top]; + col = top_index; + val = -coeff4 * att; /* Top */ + updateMatrixElement(solver_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) */ + updateMatrixElement(solver_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*/ + updateMatrixElement(solver_matrix, ptr, offset, row, col, val); + + offset = LeftStencil[StencilPosition::Center]; + col = left_index; + val = +coeff1 * arr; /* Center: (Right) -> Center: (Left) */ + updateMatrixElement(solver_matrix, ptr, offset, row, col, val); + + /* Top Right -> Bottom Left: REMOVED DUE TO ARTIFICAL 7 POINT STENCIL */ + + /* Bottom Right -> Top Left: REMOVED DUE TO ARTIFICAL 7 POINT STENCIL */ + + /* Fill matrix row of (i+1,j) */ + row = right_index; + ptr = right_nz_index; + const Stencil& RightStencil = getStencil(i_r + 1); + + offset = RightStencil[StencilPosition::Left]; + col = center_index; + val = -coeff2 * arr; /* Left */ + updateMatrixElement(solver_matrix, ptr, offset, row, col, val); + + offset = RightStencil[StencilPosition::Center]; + col = right_index; + val = +coeff2 * arr; /* Center: (Left) */ + updateMatrixElement(solver_matrix, ptr, offset, row, col, val); + + offset = RightStencil[StencilPosition::TopLeft]; + col = top_index; + val = +0.25 * art; /* Top Left */ + updateMatrixElement(solver_matrix, ptr, offset, row, col, val); + + offset = RightStencil[StencilPosition::BottomLeft]; + col = bottom_index; + val = -0.25 * art; /* Bottom Left */ + updateMatrixElement(solver_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::Top]; + col = center_index; + val = -coeff3 * att; /* Top */ + updateMatrixElement(solver_matrix, ptr, offset, row, col, val); + + offset = BottomStencil[StencilPosition::Center]; + col = bottom_index; + val = +coeff3 * att; /* Center: (Top) */ + updateMatrixElement(solver_matrix, ptr, offset, row, col, val); + + offset = BottomStencil[StencilPosition::TopRight]; + col = right_index; + val = -0.25 * art; /* Top Right */ + updateMatrixElement(solver_matrix, ptr, offset, row, col, val); + + /* TopLeft REMOVED DUE TO ARTIFICAL 7 POINT STENCIL */ + + /* Fill matrix row of (i,j+1) */ + row = top_index; + ptr = top_nz_index; + + const Stencil& TopStencil = CenterStencil; + + offset = TopStencil[StencilPosition::Bottom]; + col = center_index; + val = -coeff4 * att; /* Bottom */ + updateMatrixElement(solver_matrix, ptr, offset, row, col, val); + + offset = TopStencil[StencilPosition::Center]; + col = top_index; + val = +coeff4 * att; /* Center: (Bottom) */ + updateMatrixElement(solver_matrix, ptr, offset, row, col, val); + + offset = TopStencil[StencilPosition::BottomRight]; + col = right_index; + val = +0.25 * art; /* Bottom Right */ + updateMatrixElement(solver_matrix, ptr, offset, row, col, val); + + /* BottomLeft REMOVED DUE TO ARTIFICAL 7 POINT STENCIL */ + } + } + /* ------------------------------- */ + /* Node next to the inner boundary */ + /* ------------------------------- */ + else if (i_r == 1) { + const double h1 = grid.radialSpacing(i_r - 1); + const double h2 = grid.radialSpacing(i_r); + const double k1 = grid.angularSpacing(i_theta - 1); + const double k2 = grid.angularSpacing(i_theta); + + const double coeff1 = 0.5 * (k1 + k2) / h1; + const double coeff2 = 0.5 * (k1 + k2) / h2; + const double coeff3 = 0.5 * (h1 + h2) / k1; + const 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_nz_index = getSolverMatrixIndex(i_r, i_theta); + const int left_nz_index = getSolverMatrixIndex(i_r - 1, i_theta); + const int right_nz_index = getSolverMatrixIndex(i_r + 1, i_theta); + const int bottom_nz_index = getSolverMatrixIndex(i_r, i_theta_M1); + const int top_nz_index = getSolverMatrixIndex(i_r, i_theta_P1); + + const int center_index = grid.index(i_r, i_theta); + const int left_index = grid.index(i_r - 1, i_theta); + const int right_index = grid.index(i_r + 1, i_theta); + const int bottom_index = grid.index(i_r, i_theta_M1); + const int top_index = grid.index(i_r, i_theta_P1); + + /* Fill matrix row of (i,j) */ + row = center_index; + ptr = center_nz_index; + + const Stencil& CenterStencil = getStencil(i_r); + + offset = CenterStencil[StencilPosition::Center]; + col = center_index; + val = 0.25 * (h1 + h2) * (k1 + k2) * coeff_beta * std::fabs(detDF); /* beta_{i,j} */ + updateMatrixElement(solver_matrix, ptr, offset, row, col, val); + + /* REMOVED: Moved to the right hand side to make the matrix symmetric */ + if (!DirBC_Interior) { + offset = CenterStencil[StencilPosition::Left]; + col = left_index; + val = -coeff1 * arr; /* Left */ + updateMatrixElement(solver_matrix, ptr, offset, row, col, val); + } + + offset = CenterStencil[StencilPosition::Right]; + col = right_index; + val = -coeff2 * arr; /* Right */ + updateMatrixElement(solver_matrix, ptr, offset, row, col, val); + + offset = CenterStencil[StencilPosition::Bottom]; + col = bottom_index; + val = -coeff3 * att; /* Bottom */ + updateMatrixElement(solver_matrix, ptr, offset, row, col, val); + + offset = CenterStencil[StencilPosition::Top]; + col = top_index; + val = -coeff4 * att; /* Top */ + updateMatrixElement(solver_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) */ + updateMatrixElement(solver_matrix, ptr, offset, row, col, val); + + if (!DirBC_Interior) { /* Don't give to the inner Dirichlet boundary! */ + /* Fill matrix row of (i-1,j) */ + row = left_index; + ptr = left_nz_index; + + const Stencil& LeftStencil = getStencil(i_r - 1); + + offset = LeftStencil[StencilPosition::Right]; + col = center_index; + val = -coeff1 * arr; /* Right */ + updateMatrixElement(solver_matrix, ptr, offset, row, col, val); + + offset = LeftStencil[StencilPosition::Center]; + col = left_index; + val = +coeff1 * arr; /* Center: (Right) */ + updateMatrixElement(solver_matrix, ptr, offset, row, col, val); + + offset = LeftStencil[StencilPosition::TopRight]; + col = top_index; + val = -0.25 * art; /* Top Right */ + updateMatrixElement(solver_matrix, ptr, offset, row, col, val); + + offset = LeftStencil[StencilPosition::BottomRight]; + col = bottom_index; + val = +0.25 * art; /* Bottom Right */ + updateMatrixElement(solver_matrix, ptr, offset, row, col, val); + } + + /* Fill matrix row of (i+1,j) */ + row = right_index; + ptr = right_nz_index; + + const Stencil& RightStencil = getStencil(i_r + 1); + + offset = RightStencil[StencilPosition::Left]; + col = center_index; + val = -coeff2 * arr; /* Left */ + updateMatrixElement(solver_matrix, ptr, offset, row, col, val); + + offset = RightStencil[StencilPosition::Center]; + col = right_index; + val = +coeff2 * arr; /* Center: (Left) */ + updateMatrixElement(solver_matrix, ptr, offset, row, col, val); + + offset = RightStencil[StencilPosition::TopLeft]; + col = top_index; + val = +0.25 * art; /* Top Left */ + updateMatrixElement(solver_matrix, ptr, offset, row, col, val); + + offset = RightStencil[StencilPosition::BottomLeft]; + col = bottom_index; + val = -0.25 * art; /* Bottom Left */ + updateMatrixElement(solver_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::Top]; + col = center_index; + val = -coeff3 * att; /* Top */ + updateMatrixElement(solver_matrix, ptr, offset, row, col, val); + + offset = BottomStencil[StencilPosition::Center]; + col = bottom_index; + val = +coeff3 * att; /* Center: (Top) */ + updateMatrixElement(solver_matrix, ptr, offset, row, col, val); + + offset = BottomStencil[StencilPosition::TopRight]; + col = right_index; + val = -0.25 * art; /* Top Right */ + updateMatrixElement(solver_matrix, ptr, offset, row, col, val); + + /* REMOVED: Moved to the right hand side to make the matrix symmetric */ + if (!DirBC_Interior) { + offset = BottomStencil[StencilPosition::TopLeft]; + col = left_index; + val = +0.25 * art; /* Top Left */ + updateMatrixElement(solver_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::Bottom]; + col = center_index; + val = -coeff4 * att; /* Bottom */ + updateMatrixElement(solver_matrix, ptr, offset, row, col, val); + + offset = TopStencil[StencilPosition::Center]; + col = top_index; + val = +coeff4 * att; /* Center: (Bottom) */ + updateMatrixElement(solver_matrix, ptr, offset, row, col, val); + + offset = TopStencil[StencilPosition::BottomRight]; + col = right_index; + val = +0.25 * art; /* Bottom Right */ + updateMatrixElement(solver_matrix, ptr, offset, row, col, val); + + /* REMOVED: Moved to the right hand side to make the matrix symmetric */ + if (!DirBC_Interior) { + offset = TopStencil[StencilPosition::BottomLeft]; + col = left_index; + val = -0.25 * art; /* Bottom Left */ + updateMatrixElement(solver_matrix, ptr, offset, row, col, val); + } + } + /* ------------------------------- */ + /* Node next to the outer boundary */ + /* ------------------------------- */ + else if (i_r == grid.nr() - 2) { + const double h1 = grid.radialSpacing(i_r - 1); + const double h2 = grid.radialSpacing(i_r); + const double k1 = grid.angularSpacing(i_theta - 1); + const double k2 = grid.angularSpacing(i_theta); + + const double coeff1 = 0.5 * (k1 + k2) / h1; + const double coeff2 = 0.5 * (k1 + k2) / h2; + const double coeff3 = 0.5 * (h1 + h2) / k1; + const 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_nz_index = getSolverMatrixIndex(i_r, i_theta); + const int left_nz_index = getSolverMatrixIndex(i_r - 1, i_theta); + const int right_nz_index = getSolverMatrixIndex(i_r + 1, i_theta); + const int bottom_nz_index = getSolverMatrixIndex(i_r, i_theta_M1); + const int top_nz_index = getSolverMatrixIndex(i_r, i_theta_P1); + + const int center_index = grid.index(i_r, i_theta); + const int left_index = grid.index(i_r - 1, i_theta); + const int right_index = grid.index(i_r + 1, i_theta); + const int bottom_index = grid.index(i_r, i_theta_M1); + const int top_index = grid.index(i_r, i_theta_P1); + + /* Fill matrix row of (i,j) */ + row = center_index; + ptr = center_nz_index; + + const Stencil& CenterStencil = getStencil(i_r); + + offset = CenterStencil[StencilPosition::Center]; + col = center_index; + val = 0.25 * (h1 + h2) * (k1 + k2) * coeff_beta * std::fabs(detDF); /* beta_{i,j} */ + updateMatrixElement(solver_matrix, ptr, offset, row, col, val); + + offset = CenterStencil[StencilPosition::Left]; + col = left_index; + val = -coeff1 * arr; /* Left */ + updateMatrixElement(solver_matrix, ptr, offset, row, col, val); + + /* Right REMOVED: Moved to the right hand side to make the matrix symmetric */ + + offset = CenterStencil[StencilPosition::Bottom]; + col = bottom_index; + val = -coeff3 * att; /* Bottom */ + updateMatrixElement(solver_matrix, ptr, offset, row, col, val); + + offset = CenterStencil[StencilPosition::Top]; + col = top_index; + val = -coeff4 * att; /* Top */ + updateMatrixElement(solver_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) */ + updateMatrixElement(solver_matrix, ptr, offset, row, col, val); + + /* Fill matrix row of (i-1,j) */ + row = left_index; + ptr = left_nz_index; + + const Stencil& LeftStencil = getStencil(i_r - 1); + + offset = LeftStencil[StencilPosition::Right]; + col = center_index; + val = -coeff1 * arr; /* Right */ + updateMatrixElement(solver_matrix, ptr, offset, row, col, val); + + offset = LeftStencil[StencilPosition::Center]; + col = left_index; + val = coeff1 * arr; /* Center: (Right) */ + updateMatrixElement(solver_matrix, ptr, offset, row, col, val); + + offset = LeftStencil[StencilPosition::TopRight]; + col = top_index; + val = -0.25 * art; /* Top Right */ + updateMatrixElement(solver_matrix, ptr, offset, row, col, val); + + offset = LeftStencil[StencilPosition::BottomRight]; + col = bottom_index; + val = 0.25 * art; /* Bottom Right */ + updateMatrixElement(solver_matrix, ptr, offset, row, col, val); + + /* Fill matrix row of (i+1,j) */ + /* Don't give to the outer dirichlet boundary! */ + + /* Fill matrix row of (i,j-1) */ + row = bottom_index; + ptr = bottom_nz_index; + + const Stencil& BottomStencil = CenterStencil; + + offset = BottomStencil[StencilPosition::Top]; + col = center_index; + val = -coeff3 * att; /* Top */ + updateMatrixElement(solver_matrix, ptr, offset, row, col, val); + + offset = BottomStencil[StencilPosition::Center]; + col = bottom_index; + val = coeff3 * att; /* Center: (Top) */ + updateMatrixElement(solver_matrix, ptr, offset, row, col, val); + + /* TopRight REMOVED: Moved to the right hand side to make the matrix symmetric */ + + offset = BottomStencil[StencilPosition::TopLeft]; + col = left_index; + val = 0.25 * art; /* Top Left */ + updateMatrixElement(solver_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::Bottom]; + col = center_index; + val = -coeff4 * att; /* Bottom */ + updateMatrixElement(solver_matrix, ptr, offset, row, col, val); + + offset = TopStencil[StencilPosition::Center]; + col = top_index; + val = coeff4 * att; /* Center: (Bottom) */ + updateMatrixElement(solver_matrix, ptr, offset, row, col, val); + + /* BottomRight REMOVED: Moved to the right hand side to make the matrix symmetric */ + + offset = TopStencil[StencilPosition::BottomLeft]; + col = left_index; + val = -0.25 * art; /* Bottom Left */ + updateMatrixElement(solver_matrix, ptr, offset, row, col, val); + } + /* ------------------------------------ */ + /* Node on the outer dirichlet boundary */ + /* ------------------------------------ */ + else if (i_r == grid.nr() - 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; + + const int i_theta_M1 = grid.wrapThetaIndex(i_theta - 1); + const int i_theta_P1 = grid.wrapThetaIndex(i_theta + 1); + + const int center_nz_index = getSolverMatrixIndex(i_r, i_theta); + const int left_nz_index = getSolverMatrixIndex(i_r - 1, i_theta); + + const int center_index = grid.index(i_r, i_theta); + const int left_index = grid.index(i_r - 1, i_theta); + const int bottom_index = grid.index(i_r, i_theta_M1); + const int top_index = grid.index(i_r, i_theta_P1); + + /* Fill matrix row of (i,j) */ + row = center_index; + ptr = center_nz_index; + + const Stencil& CenterStencil = getStencil(i_r); + + offset = CenterStencil[StencilPosition::Center]; + col = center_index; + val = 1.0; + updateMatrixElement(solver_matrix, ptr, offset, row, col, val); + + /* Give value to the interior nodes! */ + /* Fill matrix row of (i-1,j) */ + row = left_index; + ptr = left_nz_index; + const Stencil& LeftStencil = getStencil(i_r - 1); + + /* Right REMOVED: Moved to the right hand side to make the matrix symmetric */ + + offset = LeftStencil[StencilPosition::Center]; + col = left_index; + val = coeff1 * arr; /* Center: (Right) */ + updateMatrixElement(solver_matrix, ptr, offset, row, col, val); + + /* TopRight REMOVED: Moved to the right hand side to make the matrix symmetric */ + + /* BottomRight REMOVED: Moved to the right hand side to make the matrix symmetric */ + } +} void DirectSolver_COO_MUMPS_Give::buildSolverMatrixCircleSection(const int i_r, SparseMatrixCOO& solver_matrix) { @@ -776,8 +783,8 @@ void DirectSolver_COO_MUMPS_Give::buildSolverMatrixCircleSection(const int i_r, level_cache_.obtainValues(i_r, i_theta, global_index, r, theta, coeff_beta, arr, att, art, detDF); // Build solver matrix at the current node - NODE_BUILD_SOLVER_MATRIX_GIVE(i_r, i_theta, r, theta, grid_, DirBC_Interior_, solver_matrix, arr, att, art, - detDF, coeff_beta); + nodeBuildSolverMatrixGive(i_r, i_theta, grid_, DirBC_Interior_, solver_matrix, arr, att, art, detDF, + coeff_beta); } } @@ -793,8 +800,8 @@ void DirectSolver_COO_MUMPS_Give::buildSolverMatrixRadialSection(const int i_the level_cache_.obtainValues(i_r, i_theta, global_index, r, theta, coeff_beta, arr, att, art, detDF); // Build solver matrix at the current node - NODE_BUILD_SOLVER_MATRIX_GIVE(i_r, i_theta, r, theta, grid_, DirBC_Interior_, solver_matrix, arr, att, art, - detDF, coeff_beta); + nodeBuildSolverMatrixGive(i_r, i_theta, grid_, DirBC_Interior_, solver_matrix, arr, att, art, detDF, + coeff_beta); } } diff --git a/src/DirectSolver/DirectSolver-CSR-LU-Give/buildSolverMatrix.cpp b/src/DirectSolver/DirectSolver-CSR-LU-Give/buildSolverMatrix.cpp index 11617de4..c192da23 100644 --- a/src/DirectSolver/DirectSolver-CSR-LU-Give/buildSolverMatrix.cpp +++ b/src/DirectSolver/DirectSolver-CSR-LU-Give/buildSolverMatrix.cpp @@ -2,733 +2,739 @@ #include "../../../include/common/geometry_helper.h" -#define UPDATE_MATRIX_ELEMENT(matrix, offset, row, col, val) \ - do { \ - matrix.row_nz_index(row, offset) = col; \ - matrix.row_nz_entry(row, offset) += val; \ - } while (0) - -#define NODE_BUILD_SOLVER_MATRIX_GIVE(i_r, i_theta, r, theta, grid, DirBC_Interior, solver_matrix, arr, att, art, \ - detDF, coeff_beta) \ - do { \ - int offset; \ - int row, col; \ - double val; \ - /* -------------------- */ \ - /* Node in the interior */ \ - /* -------------------- */ \ - if (i_r > 1 && i_r < grid.nr() - 2) { \ - const int i_theta_M1 = grid.wrapThetaIndex(i_theta - 1); \ - const int i_theta_P1 = grid.wrapThetaIndex(i_theta + 1); \ - \ - const double h1 = grid.radialSpacing(i_r - 1); \ - const double h2 = grid.radialSpacing(i_r); \ - const double k1 = grid.angularSpacing(i_theta_M1); \ - const double k2 = grid.angularSpacing(i_theta); \ - const double coeff1 = 0.5 * (k1 + k2) / h1; \ - const double coeff2 = 0.5 * (k1 + k2) / h2; \ - const double coeff3 = 0.5 * (h1 + h2) / k1; \ - const double coeff4 = 0.5 * (h1 + h2) / k2; \ - \ - const int center_index = grid.index(i_r, i_theta); \ - const int left_index = grid.index(i_r - 1, i_theta); \ - const int right_index = grid.index(i_r + 1, i_theta); \ - const int bottom_index = grid.index(i_r, i_theta_M1); \ - const int top_index = grid.index(i_r, i_theta_P1); \ - \ - /* Fill matrix row of (i,j) */ \ - row = center_index; \ - \ - const Stencil& CenterStencil = getStencil(i_r); \ - \ - offset = CenterStencil[StencilPosition::Center]; \ - col = center_index; \ - val = 0.25 * (h1 + h2) * (k1 + k2) * coeff_beta * fabs(detDF); /* beta_{i,j} */ \ - UPDATE_MATRIX_ELEMENT(solver_matrix, offset, row, col, val); \ - \ - offset = CenterStencil[StencilPosition::Left]; \ - col = left_index; \ - val = -coeff1 * arr; /* Left */ \ - UPDATE_MATRIX_ELEMENT(solver_matrix, offset, row, col, val); \ - \ - offset = CenterStencil[StencilPosition::Right]; \ - col = right_index; \ - val = -coeff2 * arr; /* Right */ \ - UPDATE_MATRIX_ELEMENT(solver_matrix, offset, row, col, val); \ - \ - offset = CenterStencil[StencilPosition::Bottom]; \ - col = bottom_index; \ - val = -coeff3 * att; /* Bottom */ \ - UPDATE_MATRIX_ELEMENT(solver_matrix, offset, row, col, val); \ - \ - offset = CenterStencil[StencilPosition::Top]; \ - col = top_index; \ - val = -coeff4 * att; /* Top */ \ - UPDATE_MATRIX_ELEMENT(solver_matrix, offset, row, col, val); \ - \ - offset = CenterStencil[StencilPosition::Center]; \ - col = center_index; \ - val = (coeff1 + coeff2) * arr + (coeff3 + coeff4) * att; /* Center: (Left, Right, Bottom, Top) */ \ - UPDATE_MATRIX_ELEMENT(solver_matrix, offset, row, col, val); \ - \ - /* Fill matrix row of (i-1,j) */ \ - row = left_index; \ - \ - const Stencil& LeftStencil = getStencil(i_r - 1); \ - \ - offset = LeftStencil[StencilPosition::Right]; \ - col = center_index; \ - val = -coeff1 * arr; /* Right */ \ - UPDATE_MATRIX_ELEMENT(solver_matrix, offset, row, col, val); \ - \ - offset = LeftStencil[StencilPosition::Center]; \ - col = left_index; \ - val = +coeff1 * arr; /* Center: (Right) */ \ - UPDATE_MATRIX_ELEMENT(solver_matrix, offset, row, col, val); \ - \ - offset = LeftStencil[StencilPosition::TopRight]; \ - col = top_index; \ - val = -0.25 * art; /* Top Right */ \ - UPDATE_MATRIX_ELEMENT(solver_matrix, offset, row, col, val); \ - \ - offset = LeftStencil[StencilPosition::BottomRight]; \ - col = bottom_index; \ - val = +0.25 * art; /* Bottom Right */ \ - UPDATE_MATRIX_ELEMENT(solver_matrix, offset, row, col, val); \ - \ - /* Fill matrix row of (i+1,j) */ \ - row = right_index; \ - \ - const Stencil& RightStencil = getStencil(i_r + 1); \ - \ - offset = RightStencil[StencilPosition::Left]; \ - col = center_index; \ - val = -coeff2 * arr; /* Left */ \ - UPDATE_MATRIX_ELEMENT(solver_matrix, offset, row, col, val); \ - \ - offset = RightStencil[StencilPosition::Center]; \ - col = right_index; \ - val = +coeff2 * arr; /* Center: (Left) */ \ - UPDATE_MATRIX_ELEMENT(solver_matrix, offset, row, col, val); \ - \ - offset = RightStencil[StencilPosition::TopLeft]; \ - col = top_index; \ - val = +0.25 * art; /* Top Left */ \ - UPDATE_MATRIX_ELEMENT(solver_matrix, offset, row, col, val); \ - \ - offset = RightStencil[StencilPosition::BottomLeft]; \ - col = bottom_index; \ - val = -0.25 * art; /* Bottom Left */ \ - UPDATE_MATRIX_ELEMENT(solver_matrix, offset, row, col, val); \ - \ - /* Fill matrix row of (i,j-1) */ \ - row = bottom_index; \ - \ - const Stencil& BottomStencil = CenterStencil; \ - \ - offset = BottomStencil[StencilPosition::Top]; \ - col = center_index; \ - val = -coeff3 * att; /* Top */ \ - UPDATE_MATRIX_ELEMENT(solver_matrix, offset, row, col, val); \ - \ - offset = BottomStencil[StencilPosition::Center]; \ - col = bottom_index; \ - val = +coeff3 * att; /* Center: (Top) */ \ - UPDATE_MATRIX_ELEMENT(solver_matrix, offset, row, col, val); \ - \ - offset = BottomStencil[StencilPosition::TopRight]; \ - col = right_index; \ - val = -0.25 * art; /* Top Right */ \ - UPDATE_MATRIX_ELEMENT(solver_matrix, offset, row, col, val); \ - \ - offset = BottomStencil[StencilPosition::TopLeft]; \ - col = left_index; \ - val = +0.25 * art; /* Top Left */ \ - UPDATE_MATRIX_ELEMENT(solver_matrix, offset, row, col, val); \ - \ - /* Fill matrix row of (i,j+1) */ \ - row = top_index; \ - \ - const Stencil& TopStencil = CenterStencil; \ - \ - offset = TopStencil[StencilPosition::Bottom]; \ - col = center_index; \ - val = -coeff4 * att; /* Bottom */ \ - UPDATE_MATRIX_ELEMENT(solver_matrix, offset, row, col, val); \ - \ - offset = TopStencil[StencilPosition::Center]; \ - col = top_index; \ - val = +coeff4 * att; /* Center: (Bottom) */ \ - UPDATE_MATRIX_ELEMENT(solver_matrix, offset, row, col, val); \ - \ - offset = TopStencil[StencilPosition::BottomRight]; \ - col = right_index; \ - val = +0.25 * art; /* Bottom Right */ \ - UPDATE_MATRIX_ELEMENT(solver_matrix, offset, row, col, val); \ - \ - offset = TopStencil[StencilPosition::BottomLeft]; \ - col = left_index; \ - val = -0.25 * art; /* Bottom Left */ \ - UPDATE_MATRIX_ELEMENT(solver_matrix, offset, row, col, val); \ - } \ - /* -------------------------- */ \ - /* Node on the inner boundary */ \ - /* -------------------------- */ \ - else if (i_r == 0) { \ - /* ------------------------------------------------ */ \ - /* Case 1: Dirichlet boundary on the inner boundary */ \ - /* ------------------------------------------------ */ \ - if (DirBC_Interior) { \ - const double h2 = grid.radialSpacing(i_r); \ - const double k1 = grid.angularSpacing(i_theta - 1); \ - const double k2 = grid.angularSpacing(i_theta); \ - const double coeff2 = 0.5 * (k1 + k2) / h2; \ - \ - const int i_theta_M1 = grid.wrapThetaIndex(i_theta - 1); \ - const int i_theta_P1 = grid.wrapThetaIndex(i_theta + 1); \ - \ - const int center_index = grid.index(i_r, i_theta); \ - const int right_index = grid.index(i_r + 1, i_theta); \ - const int bottom_index = grid.index(i_r, i_theta_M1); \ - const int top_index = grid.index(i_r, i_theta_P1); \ - \ - /* Fill matrix row of (i,j) */ \ - row = center_index; \ - \ - const Stencil& CenterStencil = getStencil(i_r); \ - \ - offset = CenterStencil[StencilPosition::Center]; \ - col = center_index; \ - val = 1.0; \ - UPDATE_MATRIX_ELEMENT(solver_matrix, offset, row, col, val); \ - \ - /* Fill matrix row of (i+1,j) */ \ - row = right_index; \ - \ - const Stencil& RightStencil = getStencil(i_r + 1); \ - \ - offset = RightStencil[StencilPosition::Left]; \ - col = center_index; \ - val = -coeff2 * arr; /* Left */ \ - UPDATE_MATRIX_ELEMENT(solver_matrix, offset, row, col, val); \ - \ - offset = RightStencil[StencilPosition::Center]; \ - col = right_index; \ - val = +coeff2 * arr; /* Center: (Left) */ \ - UPDATE_MATRIX_ELEMENT(solver_matrix, offset, row, col, val); \ - \ - offset = RightStencil[StencilPosition::TopLeft]; \ - col = top_index; \ - val = +0.25 * art; /* Top Left */ \ - UPDATE_MATRIX_ELEMENT(solver_matrix, offset, row, col, val); \ - \ - offset = RightStencil[StencilPosition::BottomLeft]; \ - col = bottom_index; \ - val = -0.25 * art; /* Bottom Left */ \ - UPDATE_MATRIX_ELEMENT(solver_matrix, offset, row, col, val); \ - } \ - 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. */ \ - const int i_theta_M1 = grid.wrapThetaIndex(i_theta - 1); \ - const int i_theta_P1 = grid.wrapThetaIndex(i_theta + 1); \ - \ - assert(grid_.ntheta() % 2 == 0); \ - const int i_theta_AcrossOrigin = grid.wrapThetaIndex(i_theta + grid.ntheta() / 2); \ - \ - double h1 = 2.0 * grid.radius(0); \ - double h2 = grid.radialSpacing(i_r); \ - double k1 = grid.angularSpacing(i_theta_M1); \ - 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 center_index = grid.index(i_r, i_theta); \ - const int left_index = grid.index(i_r, i_theta_AcrossOrigin); \ - const int right_index = grid.index(i_r + 1, i_theta); \ - const int bottom_index = grid.index(i_r, i_theta_M1); \ - const int top_index = grid.index(i_r, i_theta_P1); \ - \ - /* Fill matrix row of (i,j) */ \ - row = center_index; \ - \ - const Stencil& CenterStencil = getStencil(i_r); \ - \ - offset = CenterStencil[StencilPosition::Center]; \ - col = center_index; \ - val = 0.25 * (h1 + h2) * (k1 + k2) * coeff_beta * fabs(detDF); /* beta_{i,j} */ \ - UPDATE_MATRIX_ELEMENT(solver_matrix, offset, row, col, val); \ - \ - offset = CenterStencil[StencilPosition::Left]; \ - col = left_index; \ - val = -coeff1 * arr; /* Left */ \ - UPDATE_MATRIX_ELEMENT(solver_matrix, offset, row, col, val); \ - \ - offset = CenterStencil[StencilPosition::Right]; \ - col = right_index; \ - val = -coeff2 * arr; /* Right */ \ - UPDATE_MATRIX_ELEMENT(solver_matrix, offset, row, col, val); \ - \ - offset = CenterStencil[StencilPosition::Bottom]; \ - col = bottom_index; \ - val = -coeff3 * att; /* Bottom */ \ - UPDATE_MATRIX_ELEMENT(solver_matrix, offset, row, col, val); \ - \ - offset = CenterStencil[StencilPosition::Top]; \ - col = top_index; \ - val = -coeff4 * att; /* Top */ \ - UPDATE_MATRIX_ELEMENT(solver_matrix, offset, row, col, val); \ - \ - offset = CenterStencil[StencilPosition::Center]; \ - col = center_index; \ - val = (coeff1 + coeff2) * arr + (coeff3 + coeff4) * att; /* Center: (Left, Right, Bottom, Top) */ \ - UPDATE_MATRIX_ELEMENT(solver_matrix, 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; \ - \ - const Stencil& LeftStencil = CenterStencil; \ - \ - offset = LeftStencil[StencilPosition::Left]; \ - col = center_index; \ - val = -coeff1 * arr; /* Right -> Left*/ \ - UPDATE_MATRIX_ELEMENT(solver_matrix, offset, row, col, val); \ - \ - offset = LeftStencil[StencilPosition::Center]; \ - col = left_index; \ - val = +coeff1 * arr; /* Center: (Right) -> Center: (Left) */ \ - UPDATE_MATRIX_ELEMENT(solver_matrix, offset, row, col, val); \ - \ - /* Top Right -> Bottom Left: REMOVED DUE TO ARTIFICAL 7 POINT STENCIL */ \ - \ - /* Bottom Right -> Top Left: REMOVED DUE TO ARTIFICAL 7 POINT STENCIL */ \ - \ - /* Fill matrix row of (i+1,j) */ \ - row = right_index; \ - \ - const Stencil& RightStencil = getStencil(i_r + 1); \ - \ - offset = RightStencil[StencilPosition::Left]; \ - col = center_index; \ - val = -coeff2 * arr; /* Left */ \ - UPDATE_MATRIX_ELEMENT(solver_matrix, offset, row, col, val); \ - \ - offset = RightStencil[StencilPosition::Center]; \ - col = right_index; \ - val = +coeff2 * arr; /* Center: (Left) */ \ - UPDATE_MATRIX_ELEMENT(solver_matrix, offset, row, col, val); \ - \ - offset = RightStencil[StencilPosition::TopLeft]; \ - col = top_index; \ - val = +0.25 * art; /* Top Left */ \ - UPDATE_MATRIX_ELEMENT(solver_matrix, offset, row, col, val); \ - \ - offset = RightStencil[StencilPosition::BottomLeft]; \ - col = bottom_index; \ - val = -0.25 * art; /* Bottom Left */ \ - UPDATE_MATRIX_ELEMENT(solver_matrix, offset, row, col, val); \ - \ - /* Fill matrix row of (i,j-1) */ \ - row = bottom_index; \ - \ - const Stencil& BottomStencil = CenterStencil; \ - \ - offset = BottomStencil[StencilPosition::Top]; \ - col = center_index; \ - val = -coeff3 * att; /* Top */ \ - UPDATE_MATRIX_ELEMENT(solver_matrix, offset, row, col, val); \ - \ - offset = BottomStencil[StencilPosition::Center]; \ - col = bottom_index; \ - val = +coeff3 * att; /* Center: (Top) */ \ - UPDATE_MATRIX_ELEMENT(solver_matrix, offset, row, col, val); \ - \ - offset = BottomStencil[StencilPosition::TopRight]; \ - col = right_index; \ - val = -0.25 * art; /* Top Right */ \ - UPDATE_MATRIX_ELEMENT(solver_matrix, offset, row, col, val); \ - \ - /* TopLeft REMOVED DUE TO ARTIFICAL 7 POINT STENCIL */ \ - \ - /* Fill matrix row of (i,j+1) */ \ - row = top_index; \ - \ - const Stencil& TopStencil = CenterStencil; \ - \ - offset = TopStencil[StencilPosition::Bottom]; \ - col = center_index; \ - val = -coeff4 * att; /* Bottom */ \ - UPDATE_MATRIX_ELEMENT(solver_matrix, offset, row, col, val); \ - \ - offset = TopStencil[StencilPosition::Center]; \ - col = top_index; \ - val = +coeff4 * att; /* Center: (Bottom) */ \ - UPDATE_MATRIX_ELEMENT(solver_matrix, offset, row, col, val); \ - \ - offset = TopStencil[StencilPosition::BottomRight]; \ - col = right_index; \ - val = +0.25 * art; /* Bottom Right */ \ - UPDATE_MATRIX_ELEMENT(solver_matrix, offset, row, col, val); \ - \ - /* BottomLeft REMOVED DUE TO ARTIFICAL 7 POINT STENCIL */ \ - } \ - } \ - /* ------------------------------- */ \ - /* Node next to the inner boundary */ \ - /* ------------------------------- */ \ - else if (i_r == 1) { \ - const double h1 = grid.radialSpacing(i_r - 1); \ - const double h2 = grid.radialSpacing(i_r); \ - const double k1 = grid.angularSpacing(i_theta - 1); \ - const double k2 = grid.angularSpacing(i_theta); \ - const double coeff1 = 0.5 * (k1 + k2) / h1; \ - const double coeff2 = 0.5 * (k1 + k2) / h2; \ - const double coeff3 = 0.5 * (h1 + h2) / k1; \ - const 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 = grid.index(i_r, i_theta); \ - const int left_index = grid.index(i_r - 1, i_theta); \ - const int right_index = grid.index(i_r + 1, i_theta); \ - const int bottom_index = grid.index(i_r, i_theta_M1); \ - const int top_index = grid.index(i_r, i_theta_P1); \ - \ - /* Fill matrix row of (i,j) */ \ - row = center_index; \ - \ - const Stencil& CenterStencil = getStencil(i_r); \ - \ - offset = CenterStencil[StencilPosition::Center]; \ - col = center_index; \ - val = 0.25 * (h1 + h2) * (k1 + k2) * coeff_beta * fabs(detDF); /* beta_{i,j} */ \ - UPDATE_MATRIX_ELEMENT(solver_matrix, offset, row, col, val); \ - \ - offset = CenterStencil[StencilPosition::Left]; \ - col = left_index; \ - val = -coeff1 * arr; /* Left */ \ - UPDATE_MATRIX_ELEMENT(solver_matrix, offset, row, col, val); \ - \ - offset = CenterStencil[StencilPosition::Right]; \ - col = right_index; \ - val = -coeff2 * arr; /* Right */ \ - UPDATE_MATRIX_ELEMENT(solver_matrix, offset, row, col, val); \ - \ - offset = CenterStencil[StencilPosition::Bottom]; \ - col = bottom_index; \ - val = -coeff3 * att; /* Bottom */ \ - UPDATE_MATRIX_ELEMENT(solver_matrix, offset, row, col, val); \ - \ - offset = CenterStencil[StencilPosition::Top]; \ - col = top_index; \ - val = -coeff4 * att; /* Top */ \ - UPDATE_MATRIX_ELEMENT(solver_matrix, offset, row, col, val); \ - \ - offset = CenterStencil[StencilPosition::Center]; \ - col = center_index; \ - val = (coeff1 + coeff2) * arr + (coeff3 + coeff4) * att; /* Center: (Left, Right, Bottom, Top) */ \ - UPDATE_MATRIX_ELEMENT(solver_matrix, offset, row, col, val); \ - \ - if (!DirBC_Interior) { /* Don't give to the inner Dirichlet boundary! */ \ - /* Fill matrix row of (i-1,j) */ \ - row = left_index; \ - \ - const Stencil& LeftStencil = getStencil(i_r - 1); \ - \ - offset = LeftStencil[StencilPosition::Right]; \ - col = center_index; \ - val = -coeff1 * arr; /* Right */ \ - UPDATE_MATRIX_ELEMENT(solver_matrix, offset, row, col, val); \ - \ - offset = LeftStencil[StencilPosition::Center]; \ - col = left_index; \ - val = +coeff1 * arr; /* Center: (Right) */ \ - UPDATE_MATRIX_ELEMENT(solver_matrix, offset, row, col, val); \ - \ - offset = LeftStencil[StencilPosition::TopRight]; \ - col = top_index; \ - val = -0.25 * art; /* Top Right */ \ - UPDATE_MATRIX_ELEMENT(solver_matrix, offset, row, col, val); \ - \ - offset = LeftStencil[StencilPosition::BottomRight]; \ - col = bottom_index; \ - val = +0.25 * art; /* Bottom Right */ \ - UPDATE_MATRIX_ELEMENT(solver_matrix, offset, row, col, val); \ - } \ - \ - /* Fill matrix row of (i+1,j) */ \ - row = right_index; \ - \ - const Stencil& RightStencil = getStencil(i_r + 1); \ - \ - offset = RightStencil[StencilPosition::Left]; \ - col = center_index; \ - val = -coeff2 * arr; /* Left */ \ - UPDATE_MATRIX_ELEMENT(solver_matrix, offset, row, col, val); \ - \ - offset = RightStencil[StencilPosition::Center]; \ - col = right_index; \ - val = +coeff2 * arr; /* Center: (Left) */ \ - UPDATE_MATRIX_ELEMENT(solver_matrix, offset, row, col, val); \ - \ - offset = RightStencil[StencilPosition::TopLeft]; \ - col = top_index; \ - val = +0.25 * art; /* Top Left */ \ - UPDATE_MATRIX_ELEMENT(solver_matrix, offset, row, col, val); \ - \ - offset = RightStencil[StencilPosition::BottomLeft]; \ - col = bottom_index; \ - val = -0.25 * art; /* Bottom Left */ \ - UPDATE_MATRIX_ELEMENT(solver_matrix, offset, row, col, val); \ - \ - /* Fill matrix row of (i,j-1) */ \ - row = bottom_index; \ - \ - const Stencil& BottomStencil = CenterStencil; \ - \ - offset = BottomStencil[StencilPosition::Top]; \ - col = center_index; \ - val = -coeff3 * att; /* Top */ \ - UPDATE_MATRIX_ELEMENT(solver_matrix, offset, row, col, val); \ - \ - offset = BottomStencil[StencilPosition::Center]; \ - col = bottom_index; \ - val = +coeff3 * att; /* Center: (Top) */ \ - UPDATE_MATRIX_ELEMENT(solver_matrix, offset, row, col, val); \ - \ - offset = BottomStencil[StencilPosition::TopRight]; \ - col = right_index; \ - val = -0.25 * art; /* Top Right */ \ - UPDATE_MATRIX_ELEMENT(solver_matrix, offset, row, col, val); \ - \ - offset = BottomStencil[StencilPosition::TopLeft]; \ - col = left_index; \ - val = +0.25 * art; /* Top Left */ \ - UPDATE_MATRIX_ELEMENT(solver_matrix, offset, row, col, val); \ - \ - /* Fill matrix row of (i,j+1) */ \ - row = top_index; \ - \ - const Stencil& TopStencil = CenterStencil; \ - \ - offset = TopStencil[StencilPosition::Bottom]; \ - col = center_index; \ - val = -coeff4 * att; /* Bottom */ \ - UPDATE_MATRIX_ELEMENT(solver_matrix, offset, row, col, val); \ - \ - offset = TopStencil[StencilPosition::Center]; \ - col = top_index; \ - val = +coeff4 * att; /* Center: (Bottom) */ \ - UPDATE_MATRIX_ELEMENT(solver_matrix, offset, row, col, val); \ - \ - offset = TopStencil[StencilPosition::BottomRight]; \ - col = right_index; \ - val = +0.25 * art; /* Bottom Right */ \ - UPDATE_MATRIX_ELEMENT(solver_matrix, offset, row, col, val); \ - \ - offset = TopStencil[StencilPosition::BottomLeft]; \ - col = left_index; \ - val = -0.25 * art; /* Bottom Left */ \ - UPDATE_MATRIX_ELEMENT(solver_matrix, offset, row, col, val); \ - } \ - /* ------------------------------- */ \ - /* Node next to the outer boundary */ \ - /* ------------------------------- */ \ - else if (i_r == grid.nr() - 2) { \ - const double h1 = grid.radialSpacing(i_r - 1); \ - const double h2 = grid.radialSpacing(i_r); \ - const double k1 = grid.angularSpacing(i_theta - 1); \ - const double k2 = grid.angularSpacing(i_theta); \ - const double coeff1 = 0.5 * (k1 + k2) / h1; \ - const double coeff2 = 0.5 * (k1 + k2) / h2; \ - const double coeff3 = 0.5 * (h1 + h2) / k1; \ - const 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 = grid.index(i_r, i_theta); \ - const int left_index = grid.index(i_r - 1, i_theta); \ - const int right_index = grid.index(i_r + 1, i_theta); \ - const int bottom_index = grid.index(i_r, i_theta_M1); \ - const int top_index = grid.index(i_r, i_theta_P1); \ - \ - /* Fill matrix row of (i,j) */ \ - row = center_index; \ - \ - const Stencil& CenterStencil = getStencil(i_r); \ - \ - offset = CenterStencil[StencilPosition::Center]; \ - col = center_index; \ - val = 0.25 * (h1 + h2) * (k1 + k2) * coeff_beta * fabs(detDF); /* beta_{i,j} */ \ - UPDATE_MATRIX_ELEMENT(solver_matrix, offset, row, col, val); \ - \ - offset = CenterStencil[StencilPosition::Left]; \ - col = left_index; \ - val = -coeff1 * arr; /* Left */ \ - UPDATE_MATRIX_ELEMENT(solver_matrix, offset, row, col, val); \ - \ - offset = CenterStencil[StencilPosition::Right]; \ - col = right_index; \ - val = -coeff2 * arr; /* Right */ \ - UPDATE_MATRIX_ELEMENT(solver_matrix, offset, row, col, val); \ - \ - offset = CenterStencil[StencilPosition::Bottom]; \ - col = bottom_index; \ - val = -coeff3 * att; /* Bottom */ \ - UPDATE_MATRIX_ELEMENT(solver_matrix, offset, row, col, val); \ - \ - offset = CenterStencil[StencilPosition::Top]; \ - col = top_index; \ - val = -coeff4 * att; /* Top */ \ - UPDATE_MATRIX_ELEMENT(solver_matrix, offset, row, col, val); \ - \ - offset = CenterStencil[StencilPosition::Center]; \ - col = center_index; \ - val = (coeff1 + coeff2) * arr + (coeff3 + coeff4) * att; /* Center: (Left, Right, Bottom, Top) */ \ - UPDATE_MATRIX_ELEMENT(solver_matrix, offset, row, col, val); \ - \ - /* Fill matrix row of (i-1,j) */ \ - row = left_index; \ - \ - const Stencil& LeftStencil = getStencil(i_r - 1); \ - \ - offset = LeftStencil[StencilPosition::Right]; \ - col = center_index; \ - val = -coeff1 * arr; /* Right */ \ - UPDATE_MATRIX_ELEMENT(solver_matrix, offset, row, col, val); \ - \ - offset = LeftStencil[StencilPosition::Center]; \ - col = left_index; \ - val = coeff1 * arr; /* Center: (Right) */ \ - UPDATE_MATRIX_ELEMENT(solver_matrix, offset, row, col, val); \ - \ - offset = LeftStencil[StencilPosition::TopRight]; \ - col = top_index; \ - val = -0.25 * art; /* Top Right */ \ - UPDATE_MATRIX_ELEMENT(solver_matrix, offset, row, col, val); \ - \ - offset = LeftStencil[StencilPosition::BottomRight]; \ - col = bottom_index; \ - val = 0.25 * art; /* Bottom Right */ \ - UPDATE_MATRIX_ELEMENT(solver_matrix, offset, row, col, val); \ - \ - /* Fill matrix row of (i+1,j) */ \ - /* Don't give to the outer dirichlet boundary! */ \ - \ - /* Fill matrix row of (i,j-1) */ \ - row = bottom_index; \ - \ - const Stencil& BottomStencil = CenterStencil; \ - \ - offset = BottomStencil[StencilPosition::Top]; \ - col = center_index; \ - val = -coeff3 * att; /* Top */ \ - UPDATE_MATRIX_ELEMENT(solver_matrix, offset, row, col, val); \ - \ - offset = BottomStencil[StencilPosition::Center]; \ - col = bottom_index; \ - val = coeff3 * att; /* Center: (Top) */ \ - UPDATE_MATRIX_ELEMENT(solver_matrix, offset, row, col, val); \ - \ - offset = BottomStencil[StencilPosition::TopRight]; \ - col = right_index; \ - val = -0.25 * art; /* Top Right */ \ - UPDATE_MATRIX_ELEMENT(solver_matrix, offset, row, col, val); \ - \ - offset = BottomStencil[StencilPosition::TopLeft]; \ - col = left_index; \ - val = 0.25 * art; /* Top Left */ \ - UPDATE_MATRIX_ELEMENT(solver_matrix, offset, row, col, val); \ - \ - /* Fill matrix row of (i,j+1) */ \ - row = top_index; \ - \ - const Stencil& TopStencil = CenterStencil; \ - \ - offset = TopStencil[StencilPosition::Bottom]; \ - col = center_index; \ - val = -coeff4 * att; /* Bottom */ \ - UPDATE_MATRIX_ELEMENT(solver_matrix, offset, row, col, val); \ - \ - offset = TopStencil[StencilPosition::Center]; \ - col = top_index; \ - val = coeff4 * att; /* Center: (Bottom) */ \ - UPDATE_MATRIX_ELEMENT(solver_matrix, offset, row, col, val); \ - \ - offset = TopStencil[StencilPosition::BottomRight]; \ - col = right_index; \ - val = +0.25 * art; /* Bottom Right */ \ - UPDATE_MATRIX_ELEMENT(solver_matrix, offset, row, col, val); \ - \ - offset = TopStencil[StencilPosition::BottomLeft]; \ - col = left_index; \ - val = -0.25 * art; /* Bottom Left */ \ - UPDATE_MATRIX_ELEMENT(solver_matrix, offset, row, col, val); \ - } \ - /* ------------------------------------ */ \ - /* Node on the outer dirichlet boundary */ \ - /* ------------------------------------ */ \ - else if (i_r == grid.nr() - 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; \ - \ - const int i_theta_M1 = grid.wrapThetaIndex(i_theta - 1); \ - const int i_theta_P1 = grid.wrapThetaIndex(i_theta + 1); \ - \ - const int center_index = grid.index(i_r, i_theta); \ - const int left_index = grid.index(i_r - 1, i_theta); \ - const int bottom_index = grid.index(i_r, i_theta_M1); \ - const int top_index = grid.index(i_r, i_theta_P1); \ - \ - /* Fill matrix row of (i,j) */ \ - row = center_index; \ - \ - const Stencil& CenterStencil = getStencil(i_r); \ - \ - offset = CenterStencil[StencilPosition::Center]; \ - col = center_index; \ - val = 1.0; \ - UPDATE_MATRIX_ELEMENT(solver_matrix, offset, row, col, val); \ - \ - /* Give value to the interior nodes! */ \ - /* Fill matrix row of (i-1,j) */ \ - row = left_index; \ - \ - const Stencil& LeftStencil = getStencil(i_r - 1); \ - \ - offset = LeftStencil[StencilPosition::Right]; \ - col = center_index; \ - val = -coeff1 * arr; /* Right */ \ - UPDATE_MATRIX_ELEMENT(solver_matrix, offset, row, col, val); \ - \ - offset = LeftStencil[StencilPosition::Center]; \ - col = left_index; \ - val = +coeff1 * arr; /* Center: (Right) */ \ - UPDATE_MATRIX_ELEMENT(solver_matrix, offset, row, col, val); \ - \ - offset = LeftStencil[StencilPosition::TopRight]; \ - col = top_index; \ - val = -0.25 * art; /* Top Right */ \ - UPDATE_MATRIX_ELEMENT(solver_matrix, offset, row, col, val); \ - \ - offset = LeftStencil[StencilPosition::BottomRight]; \ - col = bottom_index; \ - val = +0.25 * art; /* Bottom Right */ \ - UPDATE_MATRIX_ELEMENT(solver_matrix, offset, row, col, val); \ - } \ - } while (0) +inline void updateMatrixElement(SparseMatrixCSR& matrix, int offset, int row, int col, double val) +{ + matrix.row_nz_index(row, offset) = col; + matrix.row_nz_entry(row, offset) += val; +} + +void DirectSolver_CSR_LU_Give::nodeBuildSolverMatrixGive(int i_r, int i_theta, const PolarGrid& grid, + const bool DirBC_Interior, + SparseMatrixCSR& solver_matrix, double arr, double att, + double art, double detDF, double coeff_beta) +{ + int offset; + int row, col; + double val; + /* -------------------- */ + /* Node in the interior */ + /* -------------------- */ + if (i_r > 1 && i_r < grid.nr() - 2) { + const int i_theta_M1 = grid.wrapThetaIndex(i_theta - 1); + const int i_theta_P1 = grid.wrapThetaIndex(i_theta + 1); + + const double h1 = grid.radialSpacing(i_r - 1); + const double h2 = grid.radialSpacing(i_r); + const double k1 = grid.angularSpacing(i_theta_M1); + const double k2 = grid.angularSpacing(i_theta); + + const double coeff1 = 0.5 * (k1 + k2) / h1; + const double coeff2 = 0.5 * (k1 + k2) / h2; + const double coeff3 = 0.5 * (h1 + h2) / k1; + const double coeff4 = 0.5 * (h1 + h2) / k2; + + const int center_index = grid.index(i_r, i_theta); + const int left_index = grid.index(i_r - 1, i_theta); + const int right_index = grid.index(i_r + 1, i_theta); + const int bottom_index = grid.index(i_r, i_theta_M1); + const int top_index = grid.index(i_r, i_theta_P1); + + /* Fill matrix row of (i,j) */ + row = center_index; + + const Stencil& CenterStencil = getStencil(i_r); + + offset = CenterStencil[StencilPosition::Center]; + col = center_index; + val = 0.25 * (h1 + h2) * (k1 + k2) * coeff_beta * std::fabs(detDF); /* beta_{i,j} */ + updateMatrixElement(solver_matrix, offset, row, col, val); + + offset = CenterStencil[StencilPosition::Left]; + col = left_index; + val = -coeff1 * arr; /* Left */ + updateMatrixElement(solver_matrix, offset, row, col, val); + + offset = CenterStencil[StencilPosition::Right]; + col = right_index; + val = -coeff2 * arr; /* Right */ + updateMatrixElement(solver_matrix, offset, row, col, val); + + offset = CenterStencil[StencilPosition::Bottom]; + col = bottom_index; + val = -coeff3 * att; /* Bottom */ + updateMatrixElement(solver_matrix, offset, row, col, val); + + offset = CenterStencil[StencilPosition::Top]; + col = top_index; + val = -coeff4 * att; /* Top */ + updateMatrixElement(solver_matrix, offset, row, col, val); + + offset = CenterStencil[StencilPosition::Center]; + col = center_index; + val = (coeff1 + coeff2) * arr + (coeff3 + coeff4) * att; /* Center: (Left, Right, Bottom, Top) */ + updateMatrixElement(solver_matrix, offset, row, col, val); + + /* Fill matrix row of (i-1,j) */ + row = left_index; + + const Stencil& LeftStencil = getStencil(i_r - 1); + + offset = LeftStencil[StencilPosition::Right]; + col = center_index; + val = -coeff1 * arr; /* Right */ + updateMatrixElement(solver_matrix, offset, row, col, val); + + offset = LeftStencil[StencilPosition::Center]; + col = left_index; + val = +coeff1 * arr; /* Center: (Right) */ + updateMatrixElement(solver_matrix, offset, row, col, val); + + offset = LeftStencil[StencilPosition::TopRight]; + col = top_index; + val = -0.25 * art; /* Top Right */ + updateMatrixElement(solver_matrix, offset, row, col, val); + + offset = LeftStencil[StencilPosition::BottomRight]; + col = bottom_index; + val = +0.25 * art; /* Bottom Right */ + updateMatrixElement(solver_matrix, offset, row, col, val); + + /* Fill matrix row of (i+1,j) */ + row = right_index; + + const Stencil& RightStencil = getStencil(i_r + 1); + + offset = RightStencil[StencilPosition::Left]; + col = center_index; + val = -coeff2 * arr; /* Left */ + updateMatrixElement(solver_matrix, offset, row, col, val); + + offset = RightStencil[StencilPosition::Center]; + col = right_index; + val = +coeff2 * arr; /* Center: (Left) */ + updateMatrixElement(solver_matrix, offset, row, col, val); + + offset = RightStencil[StencilPosition::TopLeft]; + col = top_index; + val = +0.25 * art; /* Top Left */ + updateMatrixElement(solver_matrix, offset, row, col, val); + + offset = RightStencil[StencilPosition::BottomLeft]; + col = bottom_index; + val = -0.25 * art; /* Bottom Left */ + updateMatrixElement(solver_matrix, offset, row, col, val); + + /* Fill matrix row of (i,j-1) */ + row = bottom_index; + + const Stencil& BottomStencil = CenterStencil; + + offset = BottomStencil[StencilPosition::Top]; + col = center_index; + val = -coeff3 * att; /* Top */ + updateMatrixElement(solver_matrix, offset, row, col, val); + + offset = BottomStencil[StencilPosition::Center]; + col = bottom_index; + val = +coeff3 * att; /* Center: (Top) */ + updateMatrixElement(solver_matrix, offset, row, col, val); + + offset = BottomStencil[StencilPosition::TopRight]; + col = right_index; + val = -0.25 * art; /* Top Right */ + updateMatrixElement(solver_matrix, offset, row, col, val); + + offset = BottomStencil[StencilPosition::TopLeft]; + col = left_index; + val = +0.25 * art; /* Top Left */ + updateMatrixElement(solver_matrix, offset, row, col, val); + + /* Fill matrix row of (i,j+1) */ + row = top_index; + + const Stencil& TopStencil = CenterStencil; + + offset = TopStencil[StencilPosition::Bottom]; + col = center_index; + val = -coeff4 * att; /* Bottom */ + updateMatrixElement(solver_matrix, offset, row, col, val); + + offset = TopStencil[StencilPosition::Center]; + col = top_index; + val = +coeff4 * att; /* Center: (Bottom) */ + updateMatrixElement(solver_matrix, offset, row, col, val); + + offset = TopStencil[StencilPosition::BottomRight]; + col = right_index; + val = +0.25 * art; /* Bottom Right */ + updateMatrixElement(solver_matrix, offset, row, col, val); + + offset = TopStencil[StencilPosition::BottomLeft]; + col = left_index; + val = -0.25 * art; /* Bottom Left */ + updateMatrixElement(solver_matrix, offset, row, col, val); + } + /* -------------------------- */ + /* Node on the inner boundary */ + /* -------------------------- */ + else if (i_r == 0) { + /* ------------------------------------------------ */ + /* Case 1: Dirichlet boundary on the inner boundary */ + /* ------------------------------------------------ */ + if (DirBC_Interior) { + const double h2 = grid.radialSpacing(i_r); + const double k1 = grid.angularSpacing(i_theta - 1); + const double k2 = grid.angularSpacing(i_theta); + + const double coeff2 = 0.5 * (k1 + k2) / h2; + + const int i_theta_M1 = grid.wrapThetaIndex(i_theta - 1); + const int i_theta_P1 = grid.wrapThetaIndex(i_theta + 1); + + const int center_index = grid.index(i_r, i_theta); + const int right_index = grid.index(i_r + 1, i_theta); + const int bottom_index = grid.index(i_r, i_theta_M1); + const int top_index = grid.index(i_r, i_theta_P1); + + /* Fill matrix row of (i,j) */ + row = center_index; + + const Stencil& CenterStencil = getStencil(i_r); + + offset = CenterStencil[StencilPosition::Center]; + col = center_index; + val = 1.0; + updateMatrixElement(solver_matrix, offset, row, col, val); + + /* Fill matrix row of (i+1,j) */ + row = right_index; + + const Stencil& RightStencil = getStencil(i_r + 1); + + offset = RightStencil[StencilPosition::Left]; + col = center_index; + val = -coeff2 * arr; /* Left */ + updateMatrixElement(solver_matrix, offset, row, col, val); + + offset = RightStencil[StencilPosition::Center]; + col = right_index; + val = +coeff2 * arr; /* Center: (Left) */ + updateMatrixElement(solver_matrix, offset, row, col, val); + + offset = RightStencil[StencilPosition::TopLeft]; + col = top_index; + val = +0.25 * art; /* Top Left */ + updateMatrixElement(solver_matrix, offset, row, col, val); + + offset = RightStencil[StencilPosition::BottomLeft]; + col = bottom_index; + val = -0.25 * art; /* Bottom Left */ + updateMatrixElement(solver_matrix, offset, row, col, val); + } + 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. */ + const int i_theta_M1 = grid.wrapThetaIndex(i_theta - 1); + const int i_theta_P1 = grid.wrapThetaIndex(i_theta + 1); + + assert(grid_.ntheta() % 2 == 0); + const int i_theta_AcrossOrigin = grid.wrapThetaIndex(i_theta + grid.ntheta() / 2); + + double h1 = 2.0 * grid.radius(0); + double h2 = grid.radialSpacing(i_r); + double k1 = grid.angularSpacing(i_theta_M1); + 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 center_index = grid.index(i_r, i_theta); + const int left_index = grid.index(i_r, i_theta_AcrossOrigin); + const int right_index = grid.index(i_r + 1, i_theta); + const int bottom_index = grid.index(i_r, i_theta_M1); + const int top_index = grid.index(i_r, i_theta_P1); + + /* Fill matrix row of (i,j) */ + row = center_index; + + const Stencil& CenterStencil = getStencil(i_r); + + offset = CenterStencil[StencilPosition::Center]; + col = center_index; + val = 0.25 * (h1 + h2) * (k1 + k2) * coeff_beta * std::fabs(detDF); /* beta_{i,j} */ + updateMatrixElement(solver_matrix, offset, row, col, val); + + offset = CenterStencil[StencilPosition::Left]; + col = left_index; + val = -coeff1 * arr; /* Left */ + updateMatrixElement(solver_matrix, offset, row, col, val); + + offset = CenterStencil[StencilPosition::Right]; + col = right_index; + val = -coeff2 * arr; /* Right */ + updateMatrixElement(solver_matrix, offset, row, col, val); + + offset = CenterStencil[StencilPosition::Bottom]; + col = bottom_index; + val = -coeff3 * att; /* Bottom */ + updateMatrixElement(solver_matrix, offset, row, col, val); + + offset = CenterStencil[StencilPosition::Top]; + col = top_index; + val = -coeff4 * att; /* Top */ + updateMatrixElement(solver_matrix, offset, row, col, val); + + offset = CenterStencil[StencilPosition::Center]; + col = center_index; + val = (coeff1 + coeff2) * arr + (coeff3 + coeff4) * att; /* Center: (Left, Right, Bottom, Top) */ + updateMatrixElement(solver_matrix, 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; + + const Stencil& LeftStencil = CenterStencil; + + offset = LeftStencil[StencilPosition::Left]; + col = center_index; + val = -coeff1 * arr; /* Right -> Left*/ + updateMatrixElement(solver_matrix, offset, row, col, val); + + offset = LeftStencil[StencilPosition::Center]; + col = left_index; + val = +coeff1 * arr; /* Center: (Right) -> Center: (Left) */ + updateMatrixElement(solver_matrix, offset, row, col, val); + + /* Top Right -> Bottom Left: REMOVED DUE TO ARTIFICAL 7 POINT STENCIL */ + + /* Bottom Right -> Top Left: REMOVED DUE TO ARTIFICAL 7 POINT STENCIL */ + + /* Fill matrix row of (i+1,j) */ + row = right_index; + + const Stencil& RightStencil = getStencil(i_r + 1); + + offset = RightStencil[StencilPosition::Left]; + col = center_index; + val = -coeff2 * arr; /* Left */ + updateMatrixElement(solver_matrix, offset, row, col, val); + + offset = RightStencil[StencilPosition::Center]; + col = right_index; + val = +coeff2 * arr; /* Center: (Left) */ + updateMatrixElement(solver_matrix, offset, row, col, val); + + offset = RightStencil[StencilPosition::TopLeft]; + col = top_index; + val = +0.25 * art; /* Top Left */ + updateMatrixElement(solver_matrix, offset, row, col, val); + + offset = RightStencil[StencilPosition::BottomLeft]; + col = bottom_index; + val = -0.25 * art; /* Bottom Left */ + updateMatrixElement(solver_matrix, offset, row, col, val); + + /* Fill matrix row of (i,j-1) */ + row = bottom_index; + + const Stencil& BottomStencil = CenterStencil; + + offset = BottomStencil[StencilPosition::Top]; + col = center_index; + val = -coeff3 * att; /* Top */ + updateMatrixElement(solver_matrix, offset, row, col, val); + + offset = BottomStencil[StencilPosition::Center]; + col = bottom_index; + val = +coeff3 * att; /* Center: (Top) */ + updateMatrixElement(solver_matrix, offset, row, col, val); + + offset = BottomStencil[StencilPosition::TopRight]; + col = right_index; + val = -0.25 * art; /* Top Right */ + updateMatrixElement(solver_matrix, offset, row, col, val); + + /* TopLeft REMOVED DUE TO ARTIFICAL 7 POINT STENCIL */ + + /* Fill matrix row of (i,j+1) */ + row = top_index; + + const Stencil& TopStencil = CenterStencil; + + offset = TopStencil[StencilPosition::Bottom]; + col = center_index; + val = -coeff4 * att; /* Bottom */ + updateMatrixElement(solver_matrix, offset, row, col, val); + + offset = TopStencil[StencilPosition::Center]; + col = top_index; + val = +coeff4 * att; /* Center: (Bottom) */ + updateMatrixElement(solver_matrix, offset, row, col, val); + + offset = TopStencil[StencilPosition::BottomRight]; + col = right_index; + val = +0.25 * art; /* Bottom Right */ + updateMatrixElement(solver_matrix, offset, row, col, val); + + /* BottomLeft REMOVED DUE TO ARTIFICAL 7 POINT STENCIL */ + } + } + /* ------------------------------- */ + /* Node next to the inner boundary */ + /* ------------------------------- */ + else if (i_r == 1) { + const double h1 = grid.radialSpacing(i_r - 1); + const double h2 = grid.radialSpacing(i_r); + const double k1 = grid.angularSpacing(i_theta - 1); + const double k2 = grid.angularSpacing(i_theta); + + const double coeff1 = 0.5 * (k1 + k2) / h1; + const double coeff2 = 0.5 * (k1 + k2) / h2; + const double coeff3 = 0.5 * (h1 + h2) / k1; + const 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 = grid.index(i_r, i_theta); + const int left_index = grid.index(i_r - 1, i_theta); + const int right_index = grid.index(i_r + 1, i_theta); + const int bottom_index = grid.index(i_r, i_theta_M1); + const int top_index = grid.index(i_r, i_theta_P1); + + /* Fill matrix row of (i,j) */ + row = center_index; + + const Stencil& CenterStencil = getStencil(i_r); + + offset = CenterStencil[StencilPosition::Center]; + col = center_index; + val = 0.25 * (h1 + h2) * (k1 + k2) * coeff_beta * std::fabs(detDF); /* beta_{i,j} */ + updateMatrixElement(solver_matrix, offset, row, col, val); + + offset = CenterStencil[StencilPosition::Left]; + col = left_index; + val = -coeff1 * arr; /* Left */ + updateMatrixElement(solver_matrix, offset, row, col, val); + + offset = CenterStencil[StencilPosition::Right]; + col = right_index; + val = -coeff2 * arr; /* Right */ + updateMatrixElement(solver_matrix, offset, row, col, val); + + offset = CenterStencil[StencilPosition::Bottom]; + col = bottom_index; + val = -coeff3 * att; /* Bottom */ + updateMatrixElement(solver_matrix, offset, row, col, val); + + offset = CenterStencil[StencilPosition::Top]; + col = top_index; + val = -coeff4 * att; /* Top */ + updateMatrixElement(solver_matrix, offset, row, col, val); + + offset = CenterStencil[StencilPosition::Center]; + col = center_index; + val = (coeff1 + coeff2) * arr + (coeff3 + coeff4) * att; /* Center: (Left, Right, Bottom, Top) */ + updateMatrixElement(solver_matrix, offset, row, col, val); + + if (!DirBC_Interior) { /* Don't give to the inner Dirichlet boundary! */ + /* Fill matrix row of (i-1,j) */ + row = left_index; + + const Stencil& LeftStencil = getStencil(i_r - 1); + + offset = LeftStencil[StencilPosition::Right]; + col = center_index; + val = -coeff1 * arr; /* Right */ + updateMatrixElement(solver_matrix, offset, row, col, val); + + offset = LeftStencil[StencilPosition::Center]; + col = left_index; + val = +coeff1 * arr; /* Center: (Right) */ + updateMatrixElement(solver_matrix, offset, row, col, val); + + offset = LeftStencil[StencilPosition::TopRight]; + col = top_index; + val = -0.25 * art; /* Top Right */ + updateMatrixElement(solver_matrix, offset, row, col, val); + + offset = LeftStencil[StencilPosition::BottomRight]; + col = bottom_index; + val = +0.25 * art; /* Bottom Right */ + updateMatrixElement(solver_matrix, offset, row, col, val); + } + + /* Fill matrix row of (i+1,j) */ + row = right_index; + + const Stencil& RightStencil = getStencil(i_r + 1); + + offset = RightStencil[StencilPosition::Left]; + col = center_index; + val = -coeff2 * arr; /* Left */ + updateMatrixElement(solver_matrix, offset, row, col, val); + + offset = RightStencil[StencilPosition::Center]; + col = right_index; + val = +coeff2 * arr; /* Center: (Left) */ + updateMatrixElement(solver_matrix, offset, row, col, val); + + offset = RightStencil[StencilPosition::TopLeft]; + col = top_index; + val = +0.25 * art; /* Top Left */ + updateMatrixElement(solver_matrix, offset, row, col, val); + + offset = RightStencil[StencilPosition::BottomLeft]; + col = bottom_index; + val = -0.25 * art; /* Bottom Left */ + updateMatrixElement(solver_matrix, offset, row, col, val); + + /* Fill matrix row of (i,j-1) */ + row = bottom_index; + + const Stencil& BottomStencil = CenterStencil; + + offset = BottomStencil[StencilPosition::Top]; + col = center_index; + val = -coeff3 * att; /* Top */ + updateMatrixElement(solver_matrix, offset, row, col, val); + + offset = BottomStencil[StencilPosition::Center]; + col = bottom_index; + val = +coeff3 * att; /* Center: (Top) */ + updateMatrixElement(solver_matrix, offset, row, col, val); + + offset = BottomStencil[StencilPosition::TopRight]; + col = right_index; + val = -0.25 * art; /* Top Right */ + updateMatrixElement(solver_matrix, offset, row, col, val); + + offset = BottomStencil[StencilPosition::TopLeft]; + col = left_index; + val = +0.25 * art; /* Top Left */ + updateMatrixElement(solver_matrix, offset, row, col, val); + + /* Fill matrix row of (i,j+1) */ + row = top_index; + + const Stencil& TopStencil = CenterStencil; + + offset = TopStencil[StencilPosition::Bottom]; + col = center_index; + val = -coeff4 * att; /* Bottom */ + updateMatrixElement(solver_matrix, offset, row, col, val); + + offset = TopStencil[StencilPosition::Center]; + col = top_index; + val = +coeff4 * att; /* Center: (Bottom) */ + updateMatrixElement(solver_matrix, offset, row, col, val); + + offset = TopStencil[StencilPosition::BottomRight]; + col = right_index; + val = +0.25 * art; /* Bottom Right */ + updateMatrixElement(solver_matrix, offset, row, col, val); + + offset = TopStencil[StencilPosition::BottomLeft]; + col = left_index; + val = -0.25 * art; /* Bottom Left */ + updateMatrixElement(solver_matrix, offset, row, col, val); + } + /* ------------------------------- */ + /* Node next to the outer boundary */ + /* ------------------------------- */ + else if (i_r == grid.nr() - 2) { + const double h1 = grid.radialSpacing(i_r - 1); + const double h2 = grid.radialSpacing(i_r); + const double k1 = grid.angularSpacing(i_theta - 1); + const double k2 = grid.angularSpacing(i_theta); + + const double coeff1 = 0.5 * (k1 + k2) / h1; + const double coeff2 = 0.5 * (k1 + k2) / h2; + const double coeff3 = 0.5 * (h1 + h2) / k1; + const 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 = grid.index(i_r, i_theta); + const int left_index = grid.index(i_r - 1, i_theta); + const int right_index = grid.index(i_r + 1, i_theta); + const int bottom_index = grid.index(i_r, i_theta_M1); + const int top_index = grid.index(i_r, i_theta_P1); + + /* Fill matrix row of (i,j) */ + row = center_index; + + const Stencil& CenterStencil = getStencil(i_r); + + offset = CenterStencil[StencilPosition::Center]; + col = center_index; + val = 0.25 * (h1 + h2) * (k1 + k2) * coeff_beta * std::fabs(detDF); /* beta_{i,j} */ + updateMatrixElement(solver_matrix, offset, row, col, val); + + offset = CenterStencil[StencilPosition::Left]; + col = left_index; + val = -coeff1 * arr; /* Left */ + updateMatrixElement(solver_matrix, offset, row, col, val); + + offset = CenterStencil[StencilPosition::Right]; + col = right_index; + val = -coeff2 * arr; /* Right */ + updateMatrixElement(solver_matrix, offset, row, col, val); + + offset = CenterStencil[StencilPosition::Bottom]; + col = bottom_index; + val = -coeff3 * att; /* Bottom */ + updateMatrixElement(solver_matrix, offset, row, col, val); + + offset = CenterStencil[StencilPosition::Top]; + col = top_index; + val = -coeff4 * att; /* Top */ + updateMatrixElement(solver_matrix, offset, row, col, val); + + offset = CenterStencil[StencilPosition::Center]; + col = center_index; + val = (coeff1 + coeff2) * arr + (coeff3 + coeff4) * att; /* Center: (Left, Right, Bottom, Top) */ + updateMatrixElement(solver_matrix, offset, row, col, val); + + /* Fill matrix row of (i-1,j) */ + row = left_index; + + const Stencil& LeftStencil = getStencil(i_r - 1); + + offset = LeftStencil[StencilPosition::Right]; + col = center_index; + val = -coeff1 * arr; /* Right */ + updateMatrixElement(solver_matrix, offset, row, col, val); + + offset = LeftStencil[StencilPosition::Center]; + col = left_index; + val = coeff1 * arr; /* Center: (Right) */ + updateMatrixElement(solver_matrix, offset, row, col, val); + + offset = LeftStencil[StencilPosition::TopRight]; + col = top_index; + val = -0.25 * art; /* Top Right */ + updateMatrixElement(solver_matrix, offset, row, col, val); + + offset = LeftStencil[StencilPosition::BottomRight]; + col = bottom_index; + val = 0.25 * art; /* Bottom Right */ + updateMatrixElement(solver_matrix, offset, row, col, val); + + /* Fill matrix row of (i+1,j) */ + /* Don't give to the outer dirichlet boundary! */ + + /* Fill matrix row of (i,j-1) */ + row = bottom_index; + + const Stencil& BottomStencil = CenterStencil; + + offset = BottomStencil[StencilPosition::Top]; + col = center_index; + val = -coeff3 * att; /* Top */ + updateMatrixElement(solver_matrix, offset, row, col, val); + + offset = BottomStencil[StencilPosition::Center]; + col = bottom_index; + val = coeff3 * att; /* Center: (Top) */ + updateMatrixElement(solver_matrix, offset, row, col, val); + + offset = BottomStencil[StencilPosition::TopRight]; + col = right_index; + val = -0.25 * art; /* Top Right */ + updateMatrixElement(solver_matrix, offset, row, col, val); + + offset = BottomStencil[StencilPosition::TopLeft]; + col = left_index; + val = 0.25 * art; /* Top Left */ + updateMatrixElement(solver_matrix, offset, row, col, val); + + /* Fill matrix row of (i,j+1) */ + row = top_index; + + const Stencil& TopStencil = CenterStencil; + + offset = TopStencil[StencilPosition::Bottom]; + col = center_index; + val = -coeff4 * att; /* Bottom */ + updateMatrixElement(solver_matrix, offset, row, col, val); + + offset = TopStencil[StencilPosition::Center]; + col = top_index; + val = coeff4 * att; /* Center: (Bottom) */ + updateMatrixElement(solver_matrix, offset, row, col, val); + + offset = TopStencil[StencilPosition::BottomRight]; + col = right_index; + val = +0.25 * art; /* Bottom Right */ + updateMatrixElement(solver_matrix, offset, row, col, val); + + offset = TopStencil[StencilPosition::BottomLeft]; + col = left_index; + val = -0.25 * art; /* Bottom Left */ + updateMatrixElement(solver_matrix, offset, row, col, val); + } + /* ------------------------------------ */ + /* Node on the outer dirichlet boundary */ + /* ------------------------------------ */ + else if (i_r == grid.nr() - 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; + + const int i_theta_M1 = grid.wrapThetaIndex(i_theta - 1); + const int i_theta_P1 = grid.wrapThetaIndex(i_theta + 1); + + const int center_index = grid.index(i_r, i_theta); + const int left_index = grid.index(i_r - 1, i_theta); + const int bottom_index = grid.index(i_r, i_theta_M1); + const int top_index = grid.index(i_r, i_theta_P1); + + /* Fill matrix row of (i,j) */ + row = center_index; + + const Stencil& CenterStencil = getStencil(i_r); + + offset = CenterStencil[StencilPosition::Center]; + col = center_index; + val = 1.0; + updateMatrixElement(solver_matrix, offset, row, col, val); + + /* Give value to the interior nodes! */ + /* Fill matrix row of (i-1,j) */ + row = left_index; + + const Stencil& LeftStencil = getStencil(i_r - 1); + + offset = LeftStencil[StencilPosition::Right]; + col = center_index; + val = -coeff1 * arr; /* Right */ + updateMatrixElement(solver_matrix, offset, row, col, val); + + offset = LeftStencil[StencilPosition::Center]; + col = left_index; + val = +coeff1 * arr; /* Center: (Right) */ + updateMatrixElement(solver_matrix, offset, row, col, val); + + offset = LeftStencil[StencilPosition::TopRight]; + col = top_index; + val = -0.25 * art; /* Top Right */ + updateMatrixElement(solver_matrix, offset, row, col, val); + + offset = LeftStencil[StencilPosition::BottomRight]; + col = bottom_index; + val = +0.25 * art; /* Bottom Right */ + updateMatrixElement(solver_matrix, offset, row, col, val); + } +} void DirectSolver_CSR_LU_Give::buildSolverMatrixCircleSection(const int i_r, SparseMatrixCSR& solver_matrix) { @@ -741,8 +747,8 @@ void DirectSolver_CSR_LU_Give::buildSolverMatrixCircleSection(const int i_r, Spa level_cache_.obtainValues(i_r, i_theta, global_index, r, theta, coeff_beta, arr, att, art, detDF); // Build solver matrix at the current node - NODE_BUILD_SOLVER_MATRIX_GIVE(i_r, i_theta, r, theta, grid_, DirBC_Interior_, solver_matrix, arr, att, art, - detDF, coeff_beta); + nodeBuildSolverMatrixGive(i_r, i_theta, grid_, DirBC_Interior_, solver_matrix, arr, att, art, detDF, + coeff_beta); } } @@ -757,8 +763,8 @@ void DirectSolver_CSR_LU_Give::buildSolverMatrixRadialSection(const int i_theta, level_cache_.obtainValues(i_r, i_theta, global_index, r, theta, coeff_beta, arr, att, art, detDF); // Build solver matrix at the current node - NODE_BUILD_SOLVER_MATRIX_GIVE(i_r, i_theta, r, theta, grid_, DirBC_Interior_, solver_matrix, arr, att, art, - detDF, coeff_beta); + nodeBuildSolverMatrixGive(i_r, i_theta, grid_, DirBC_Interior_, solver_matrix, arr, att, art, detDF, + coeff_beta); } }