diff --git a/include/DirectSolver/DirectSolver-COO-MUMPS-Take/directSolverTake.h b/include/DirectSolver/DirectSolver-COO-MUMPS-Take/directSolverTake.h index f497d19d..176adc22 100644 --- a/include/DirectSolver/DirectSolver-COO-MUMPS-Take/directSolverTake.h +++ b/include/DirectSolver/DirectSolver-COO-MUMPS-Take/directSolverTake.h @@ -86,6 +86,11 @@ class DirectSolver_COO_MUMPS_Take : public DirectSolver // Retrieves the stencil for the solver matrix at the given radial index. const Stencil& getStencil(int i_r) const; + + void nodeBuildSolverMatrixTake(int i_r, int i_theta, const PolarGrid& grid, bool DirBC_Interior, + SparseMatrixCOO& solver_matrix, ConstVector& arr, + ConstVector& att, ConstVector& art, ConstVector& detDF, + ConstVector& coeff_beta); }; #endif \ No newline at end of file diff --git a/include/DirectSolver/DirectSolver-CSR-LU-Take/directSolverTakeCustomLU.h b/include/DirectSolver/DirectSolver-CSR-LU-Take/directSolverTakeCustomLU.h index b194323f..fb89a094 100644 --- a/include/DirectSolver/DirectSolver-CSR-LU-Take/directSolverTakeCustomLU.h +++ b/include/DirectSolver/DirectSolver-CSR-LU-Take/directSolverTakeCustomLU.h @@ -57,4 +57,9 @@ class DirectSolver_CSR_LU_Take : public DirectSolver const Stencil& getStencil(int i_r) const; int getStencilSize(int global_index) const; + + void nodeBuildSolverMatrixTake(int i_r, int i_theta, const PolarGrid& grid, bool DirBC_Interior, + SparseMatrixCSR& solver_matrix, ConstVector& arr, + ConstVector& att, ConstVector& art, ConstVector& detDF, + ConstVector& coeff_beta); }; diff --git a/src/DirectSolver/DirectSolver-COO-MUMPS-Take/applySymmetryShift.cpp b/src/DirectSolver/DirectSolver-COO-MUMPS-Take/applySymmetryShift.cpp index b8b42721..a2f648f8 100644 --- a/src/DirectSolver/DirectSolver-COO-MUMPS-Take/applySymmetryShift.cpp +++ b/src/DirectSolver/DirectSolver-COO-MUMPS-Take/applySymmetryShift.cpp @@ -15,9 +15,9 @@ void DirectSolver_COO_MUMPS_Take::applySymmetryShiftInnerBoundary(Vector assert(level_cache_.cacheDensityProfileCoefficients()); assert(level_cache_.cacheDomainGeometry()); - const auto& arr = level_cache_.arr(); - const auto& att = level_cache_.att(); - const auto& art = level_cache_.art(); + ConstVector arr = level_cache_.arr(); + ConstVector att = level_cache_.att(); + ConstVector art = level_cache_.art(); const int i_r = 1; @@ -50,9 +50,9 @@ void DirectSolver_COO_MUMPS_Take::applySymmetryShiftOuterBoundary(Vector assert(level_cache_.cacheDensityProfileCoefficients()); assert(level_cache_.cacheDomainGeometry()); - const auto& arr = level_cache_.arr(); - const auto& att = level_cache_.att(); - const auto& art = level_cache_.art(); + ConstVector arr = level_cache_.arr(); + ConstVector att = level_cache_.att(); + ConstVector art = level_cache_.art(); const int i_r = grid_.nr() - 2; @@ -85,27 +85,18 @@ void DirectSolver_COO_MUMPS_Take::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 section - { - if (DirBC_Interior_) { - applySymmetryShiftInnerBoundary(x); - } + { + if (DirBC_Interior_) { + applySymmetryShiftInnerBoundary(x); } + } #pragma omp section - { - applySymmetryShiftOuterBoundary(x); - } + { + applySymmetryShiftOuterBoundary(x); } } } diff --git a/src/DirectSolver/DirectSolver-COO-MUMPS-Take/buildSolverMatrix.cpp b/src/DirectSolver/DirectSolver-COO-MUMPS-Take/buildSolverMatrix.cpp index 8a0a80e3..e58512c7 100644 --- a/src/DirectSolver/DirectSolver-COO-MUMPS-Take/buildSolverMatrix.cpp +++ b/src/DirectSolver/DirectSolver-COO-MUMPS-Take/buildSolverMatrix.cpp @@ -2,466 +2,470 @@ #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_TAKE(i_r, i_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 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); \ - const int bottom_left_index = grid.index(i_r - 1, i_theta_M1); \ - const int bottom_right_index = grid.index(i_r + 1, i_theta_M1); \ - const int top_left_index = grid.index(i_r - 1, i_theta_P1); \ - const int top_right_index = grid.index(i_r + 1, i_theta_P1); \ - \ - const double left_value = -coeff1 * (arr(center_index) + arr(left_index)); /* Left */ \ - const double right_value = -coeff2 * (arr(center_index) + arr(right_index)); /* Right */ \ - const double bottom_value = -coeff3 * (att(center_index) + att(bottom_index)); /* Bottom */ \ - const double top_value = -coeff4 * (att(center_index) + att(top_index)); /* Top */ \ - \ - const double center_value = (+0.25 * (h1 + h2) * (k1 + k2) * coeff_beta[center_index] * \ - fabs(detDF(center_index)) /* beta_{i,j} */ \ - - left_value /* Center: (Left) */ \ - - right_value /* Center: (Right) */ \ - - bottom_value /* Center: (Bottom) */ \ - - top_value /* Center: (Top) */ \ - ); \ - \ - const double bottom_left_value = -0.25 * (art(left_index) + art(bottom_index)); /* Bottom Left */ \ - const double bottom_right_value = +0.25 * (art(right_index) + art(bottom_index)); /* Bottom Right */ \ - const double top_left_value = +0.25 * (art(left_index) + art(top_index)); /* Top Left */ \ - const double top_right_value = -0.25 * (art(right_index) + art(top_index)); /* Top Right */ \ - \ - /* 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 = center_value; \ - UPDATE_MATRIX_ELEMENT(solver_matrix, ptr, offset, row, col, val); \ - \ - offset = CenterStencil[StencilPosition::Left]; \ - col = left_index; \ - val = left_value; \ - UPDATE_MATRIX_ELEMENT(solver_matrix, ptr, offset, row, col, val); \ - \ - offset = CenterStencil[StencilPosition::Right]; \ - col = right_index; \ - val = right_value; \ - UPDATE_MATRIX_ELEMENT(solver_matrix, ptr, offset, row, col, val); \ - \ - offset = CenterStencil[StencilPosition::Bottom]; \ - col = bottom_index; \ - val = bottom_value; \ - UPDATE_MATRIX_ELEMENT(solver_matrix, ptr, offset, row, col, val); \ - \ - offset = CenterStencil[StencilPosition::Top]; \ - col = top_index; \ - val = top_value; \ - UPDATE_MATRIX_ELEMENT(solver_matrix, ptr, offset, row, col, val); \ - \ - offset = CenterStencil[StencilPosition::BottomLeft]; \ - col = bottom_left_index; \ - val = bottom_left_value; \ - UPDATE_MATRIX_ELEMENT(solver_matrix, ptr, offset, row, col, val); \ - \ - offset = CenterStencil[StencilPosition::BottomRight]; \ - col = bottom_right_index; \ - val = bottom_right_value; \ - UPDATE_MATRIX_ELEMENT(solver_matrix, ptr, offset, row, col, val); \ - \ - offset = CenterStencil[StencilPosition::TopLeft]; \ - col = top_left_index; \ - val = top_left_value; \ - UPDATE_MATRIX_ELEMENT(solver_matrix, ptr, offset, row, col, val); \ - \ - offset = CenterStencil[StencilPosition::TopRight]; \ - col = top_right_index; \ - val = top_right_value; \ - 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 int center_nz_index = getSolverMatrixIndex(i_r, i_theta); \ - \ - const int center_index = grid.index(i_r, i_theta); \ - \ - /* 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); \ - } \ - 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); \ - \ - 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 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); \ - const int bottom_right_index = grid.index(i_r + 1, i_theta_M1); \ - const int top_right_index = grid.index(i_r + 1, i_theta_P1); \ - \ - const double left_value = -coeff1 * (arr(center_index) + arr(left_index)); /* Left */ \ - const double right_value = -coeff2 * (arr(center_index) + arr(right_index)); /* Right */ \ - const double bottom_value = -coeff3 * (att(center_index) + att(bottom_index)); /* Bottom */ \ - const double top_value = -coeff4 * (att(center_index) + att(top_index)); /* Top */ \ - \ - const double center_value = (+0.25 * (h1 + h2) * (k1 + k2) * coeff_beta[center_index] * \ - fabs(detDF(center_index)) /* beta_{i,j} */ \ - - left_value /* Center: (Left) */ \ - - right_value /* Center: (Right) */ \ - - bottom_value /* Center: (Bottom) */ \ - - top_value /* Center: (Top) */ \ - ); \ - \ - const double bottom_right_value = \ - +0.25 * (art(right_index) + art(bottom_index)); /* Bottom Right */ \ - const double top_right_value = -0.25 * (art(right_index) + art(top_index)); /* Top Right */ \ - \ - /* 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 = center_value; \ - UPDATE_MATRIX_ELEMENT(solver_matrix, ptr, offset, row, col, val); \ - \ - offset = CenterStencil[StencilPosition::Left]; \ - col = left_index; \ - val = left_value; \ - UPDATE_MATRIX_ELEMENT(solver_matrix, ptr, offset, row, col, val); \ - \ - offset = CenterStencil[StencilPosition::Right]; \ - col = right_index; \ - val = right_value; \ - UPDATE_MATRIX_ELEMENT(solver_matrix, ptr, offset, row, col, val); \ - \ - offset = CenterStencil[StencilPosition::Bottom]; \ - col = bottom_index; \ - val = bottom_value; \ - UPDATE_MATRIX_ELEMENT(solver_matrix, ptr, offset, row, col, val); \ - \ - offset = CenterStencil[StencilPosition::Top]; \ - col = top_index; \ - val = top_value; \ - UPDATE_MATRIX_ELEMENT(solver_matrix, ptr, offset, row, col, val); \ - \ - /* BottomLeft: REMOVED DUE TO ARTIFICAL 7 POINT STENCIL */ \ - \ - offset = CenterStencil[StencilPosition::BottomRight]; \ - col = bottom_right_index; \ - val = bottom_right_value; \ - UPDATE_MATRIX_ELEMENT(solver_matrix, ptr, offset, row, col, val); \ - \ - /* TopLeft: REMOVED DUE TO ARTIFICAL 7 POINT STENCIL */ \ - \ - offset = CenterStencil[StencilPosition::TopRight]; \ - col = top_right_index; \ - val = top_right_value; \ - UPDATE_MATRIX_ELEMENT(solver_matrix, ptr, offset, row, col, val); \ - } \ - } \ - /* ------------------------------- */ \ - /* Node next to the inner boundary */ \ - /* ------------------------------- */ \ - else if (i_r == 1) { \ - 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 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); \ - const int bottom_left_index = grid.index(i_r - 1, i_theta_M1); \ - const int bottom_right_index = grid.index(i_r + 1, i_theta_M1); \ - const int top_left_index = grid.index(i_r - 1, i_theta_P1); \ - const int top_right_index = grid.index(i_r + 1, i_theta_P1); \ - \ - const double left_value = -coeff1 * (arr(center_index) + arr(left_index)); /* Left */ \ - const double right_value = -coeff2 * (arr(center_index) + arr(right_index)); /* Right */ \ - const double bottom_value = -coeff3 * (att(center_index) + att(bottom_index)); /* Bottom */ \ - const double top_value = -coeff4 * (att(center_index) + att(top_index)); /* Top */ \ - \ - const double center_value = (+0.25 * (h1 + h2) * (k1 + k2) * coeff_beta[center_index] * \ - fabs(detDF(center_index)) /* beta_{i,j} */ \ - - left_value /* Center: (Left) */ \ - - right_value /* Center: (Right) */ \ - - bottom_value /* Center: (Bottom) */ \ - - top_value /* Center: (Top) */ \ - ); \ - \ - const double bottom_left_value = -0.25 * (art(left_index) + art(bottom_index)); /* Bottom Left */ \ - const double bottom_right_value = +0.25 * (art(right_index) + art(bottom_index)); /* Bottom Right */ \ - const double top_left_value = +0.25 * (art(left_index) + art(top_index)); /* Top Left */ \ - const double top_right_value = -0.25 * (art(right_index) + art(top_index)); /* Top Right */ \ - \ - /* 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 = center_value; \ - 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 = left_value; \ - UPDATE_MATRIX_ELEMENT(solver_matrix, ptr, offset, row, col, val); \ - } \ - \ - offset = CenterStencil[StencilPosition::Right]; \ - col = right_index; \ - val = right_value; \ - UPDATE_MATRIX_ELEMENT(solver_matrix, ptr, offset, row, col, val); \ - \ - offset = CenterStencil[StencilPosition::Bottom]; \ - col = bottom_index; \ - val = bottom_value; \ - UPDATE_MATRIX_ELEMENT(solver_matrix, ptr, offset, row, col, val); \ - \ - offset = CenterStencil[StencilPosition::Top]; \ - col = top_index; \ - val = top_value; \ - 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::BottomLeft]; \ - col = bottom_left_index; \ - val = bottom_left_value; \ - UPDATE_MATRIX_ELEMENT(solver_matrix, ptr, offset, row, col, val); \ - } \ - \ - offset = CenterStencil[StencilPosition::BottomRight]; \ - col = bottom_right_index; \ - val = bottom_right_value; \ - 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::TopLeft]; \ - col = top_left_index; \ - val = top_left_value; \ - UPDATE_MATRIX_ELEMENT(solver_matrix, ptr, offset, row, col, val); \ - } \ - \ - offset = CenterStencil[StencilPosition::TopRight]; \ - col = top_right_index; \ - val = top_right_value; \ - UPDATE_MATRIX_ELEMENT(solver_matrix, ptr, offset, row, col, val); \ - } \ - /* ------------------------------- */ \ - /* Node next to the outer boundary */ \ - /* ------------------------------- */ \ - else if (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 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); \ - const int bottom_left_index = grid.index(i_r - 1, i_theta_M1); \ - const int bottom_right_index = grid.index(i_r + 1, i_theta_M1); \ - const int top_left_index = grid.index(i_r - 1, i_theta_P1); \ - const int top_right_index = grid.index(i_r + 1, i_theta_P1); \ - \ - const double left_value = -coeff1 * (arr(center_index) + arr(left_index)); /* Left */ \ - const double right_value = -coeff2 * (arr(center_index) + arr(right_index)); /* Right */ \ - const double bottom_value = -coeff3 * (att(center_index) + att(bottom_index)); /* Bottom */ \ - const double top_value = -coeff4 * (att(center_index) + att(top_index)); /* Top */ \ - \ - const double center_value = (+0.25 * (h1 + h2) * (k1 + k2) * coeff_beta[center_index] * \ - fabs(detDF(center_index)) /* beta_{i,j} */ \ - - left_value /* Center: (Left) */ \ - - right_value /* Center: (Right) */ \ - - bottom_value /* Center: (Bottom) */ \ - - top_value /* Center: (Top) */ \ - ); \ - \ - const double bottom_left_value = -0.25 * (art(left_index) + art(bottom_index)); /* Bottom Left */ \ - const double bottom_right_value = +0.25 * (art(right_index) + art(bottom_index)); /* Bottom Right */ \ - const double top_left_value = +0.25 * (art(left_index) + art(top_index)); /* Top Left */ \ - const double top_right_value = -0.25 * (art(right_index) + art(top_index)); /* Top Right */ \ - \ - /* 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 = center_value; \ - UPDATE_MATRIX_ELEMENT(solver_matrix, ptr, offset, row, col, val); \ - \ - offset = CenterStencil[StencilPosition::Left]; \ - col = left_index; \ - val = left_value; \ - 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 = bottom_value; \ - UPDATE_MATRIX_ELEMENT(solver_matrix, ptr, offset, row, col, val); \ - \ - offset = CenterStencil[StencilPosition::Top]; \ - col = top_index; \ - val = top_value; \ - UPDATE_MATRIX_ELEMENT(solver_matrix, ptr, offset, row, col, val); \ - \ - offset = CenterStencil[StencilPosition::BottomLeft]; \ - col = bottom_left_index; \ - val = bottom_left_value; \ - UPDATE_MATRIX_ELEMENT(solver_matrix, ptr, offset, row, col, val); \ - \ - /* BottomRight REMOVED: Moved to the right hand side to make the matrix symmetric */ \ - \ - offset = CenterStencil[StencilPosition::TopLeft]; \ - col = top_left_index; \ - val = top_left_value; \ - UPDATE_MATRIX_ELEMENT(solver_matrix, ptr, offset, row, col, val); \ - \ - /* TopRight REMOVED: Moved to the right hand side to make the matrix symmetric */ \ - } \ - /* ------------------------------------ */ \ - /* Node on the outer dirichlet boundary */ \ - /* ------------------------------------ */ \ - else if (i_r == grid.nr() - 1) { \ - const int center_nz_index = getSolverMatrixIndex(i_r, i_theta); \ - \ - const int center_index = grid.index(i_r, i_theta); \ - \ - /* 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); \ - } \ - } 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_Take::nodeBuildSolverMatrixTake(int i_r, int i_theta, const PolarGrid& grid, + bool DirBC_Interior, SparseMatrixCOO& solver_matrix, + ConstVector& arr, ConstVector& att, + ConstVector& art, ConstVector& detDF, + ConstVector& coeff_beta) +{ + int ptr, offset; + int row, col; + double val; + /* -------------------- */ + /* Node in the interior */ + /* -------------------- */ + if (i_r > 1 && i_r < grid.nr() - 2) { + int i_theta_M1 = grid.wrapThetaIndex(i_theta - 1); + int i_theta_P1 = grid.wrapThetaIndex(i_theta + 1); + + double h1 = grid.radialSpacing(i_r - 1); + 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; + + int center_nz_index = getSolverMatrixIndex(i_r, i_theta); + + int center_index = grid.index(i_r, i_theta); + int left_index = grid.index(i_r - 1, i_theta); + int right_index = grid.index(i_r + 1, i_theta); + int bottom_index = grid.index(i_r, i_theta_M1); + int top_index = grid.index(i_r, i_theta_P1); + int bottom_left_index = grid.index(i_r - 1, i_theta_M1); + int bottom_right_index = grid.index(i_r + 1, i_theta_M1); + int top_left_index = grid.index(i_r - 1, i_theta_P1); + int top_right_index = grid.index(i_r + 1, i_theta_P1); + + double left_value = -coeff1 * (arr(center_index) + arr(left_index)); /* Left */ + double right_value = -coeff2 * (arr(center_index) + arr(right_index)); /* Right */ + double bottom_value = -coeff3 * (att(center_index) + att(bottom_index)); /* Bottom */ + double top_value = -coeff4 * (att(center_index) + att(top_index)); /* Top */ + + double center_value = + (+0.25 * (h1 + h2) * (k1 + k2) * coeff_beta[center_index] * fabs(detDF(center_index)) /* beta_{i,j} */ + - left_value /* Center: (Left) */ + - right_value /* Center: (Right) */ + - bottom_value /* Center: (Bottom) */ + - top_value /* Center: (Top) */ + ); + + double bottom_left_value = -0.25 * (art(left_index) + art(bottom_index)); /* Bottom Left */ + double bottom_right_value = +0.25 * (art(right_index) + art(bottom_index)); /* Bottom Right */ + double top_left_value = +0.25 * (art(left_index) + art(top_index)); /* Top Left */ + double top_right_value = -0.25 * (art(right_index) + art(top_index)); /* Top Right */ + + /* 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 = center_value; + updateMatrixElement(solver_matrix, ptr, offset, row, col, val); + + offset = CenterStencil[StencilPosition::Left]; + col = left_index; + val = left_value; + updateMatrixElement(solver_matrix, ptr, offset, row, col, val); + + offset = CenterStencil[StencilPosition::Right]; + col = right_index; + val = right_value; + updateMatrixElement(solver_matrix, ptr, offset, row, col, val); + + offset = CenterStencil[StencilPosition::Bottom]; + col = bottom_index; + val = bottom_value; + updateMatrixElement(solver_matrix, ptr, offset, row, col, val); + + offset = CenterStencil[StencilPosition::Top]; + col = top_index; + val = top_value; + updateMatrixElement(solver_matrix, ptr, offset, row, col, val); + + offset = CenterStencil[StencilPosition::BottomLeft]; + col = bottom_left_index; + val = bottom_left_value; + updateMatrixElement(solver_matrix, ptr, offset, row, col, val); + + offset = CenterStencil[StencilPosition::BottomRight]; + col = bottom_right_index; + val = bottom_right_value; + updateMatrixElement(solver_matrix, ptr, offset, row, col, val); + + offset = CenterStencil[StencilPosition::TopLeft]; + col = top_left_index; + val = top_left_value; + updateMatrixElement(solver_matrix, ptr, offset, row, col, val); + + offset = CenterStencil[StencilPosition::TopRight]; + col = top_right_index; + val = top_right_value; + 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) { + int center_nz_index = getSolverMatrixIndex(i_r, i_theta); + + int center_index = grid.index(i_r, i_theta); + + /* 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); + } + 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; + + int center_nz_index = getSolverMatrixIndex(i_r, i_theta); + + int center_index = grid.index(i_r, i_theta); + int left_index = grid.index(i_r, i_theta_AcrossOrigin); + int right_index = grid.index(i_r + 1, i_theta); + int bottom_index = grid.index(i_r, i_theta_M1); + int top_index = grid.index(i_r, i_theta_P1); + int bottom_right_index = grid.index(i_r + 1, i_theta_M1); + int top_right_index = grid.index(i_r + 1, i_theta_P1); + + double left_value = -coeff1 * (arr(center_index) + arr(left_index)); /* Left */ + double right_value = -coeff2 * (arr(center_index) + arr(right_index)); /* Right */ + double bottom_value = -coeff3 * (att(center_index) + att(bottom_index)); /* Bottom */ + double top_value = -coeff4 * (att(center_index) + att(top_index)); /* Top */ + + double center_value = + (+0.25 * (h1 + h2) * (k1 + k2) * coeff_beta[center_index] * fabs(detDF(center_index)) /* beta_{i,j} */ + - left_value /* Center: (Left) */ + - right_value /* Center: (Right) */ + - bottom_value /* Center: (Bottom) */ + - top_value /* Center: (Top) */ + ); + + double bottom_right_value = +0.25 * (art(right_index) + art(bottom_index)); /* Bottom Right */ + double top_right_value = -0.25 * (art(right_index) + art(top_index)); /* Top Right */ + + /* 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 = center_value; + updateMatrixElement(solver_matrix, ptr, offset, row, col, val); + + offset = CenterStencil[StencilPosition::Left]; + col = left_index; + val = left_value; + updateMatrixElement(solver_matrix, ptr, offset, row, col, val); + + offset = CenterStencil[StencilPosition::Right]; + col = right_index; + val = right_value; + updateMatrixElement(solver_matrix, ptr, offset, row, col, val); + + offset = CenterStencil[StencilPosition::Bottom]; + col = bottom_index; + val = bottom_value; + updateMatrixElement(solver_matrix, ptr, offset, row, col, val); + + offset = CenterStencil[StencilPosition::Top]; + col = top_index; + val = top_value; + updateMatrixElement(solver_matrix, ptr, offset, row, col, val); + + /* BottomLeft: REMOVED DUE TO ARTIFICAL 7 POINT STENCIL */ + + offset = CenterStencil[StencilPosition::BottomRight]; + col = bottom_right_index; + val = bottom_right_value; + updateMatrixElement(solver_matrix, ptr, offset, row, col, val); + + /* TopLeft: REMOVED DUE TO ARTIFICAL 7 POINT STENCIL */ + + offset = CenterStencil[StencilPosition::TopRight]; + col = top_right_index; + val = top_right_value; + updateMatrixElement(solver_matrix, ptr, offset, row, col, val); + } + } + /* ------------------------------- */ + /* Node next to the inner boundary */ + /* ------------------------------- */ + else if (i_r == 1) { + int i_theta_M1 = grid.wrapThetaIndex(i_theta - 1); + int i_theta_P1 = grid.wrapThetaIndex(i_theta + 1); + + double h1 = grid.radialSpacing(i_r - 1); + 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; + + int center_nz_index = getSolverMatrixIndex(i_r, i_theta); + + int center_index = grid.index(i_r, i_theta); + int left_index = grid.index(i_r - 1, i_theta); + int right_index = grid.index(i_r + 1, i_theta); + int bottom_index = grid.index(i_r, i_theta_M1); + int top_index = grid.index(i_r, i_theta_P1); + int bottom_left_index = grid.index(i_r - 1, i_theta_M1); + const int bottom_right_index = grid.index(i_r + 1, i_theta_M1); + int top_left_index = grid.index(i_r - 1, i_theta_P1); + int top_right_index = grid.index(i_r + 1, i_theta_P1); + + double left_value = -coeff1 * (arr(center_index) + arr(left_index)); /* Left */ + double right_value = -coeff2 * (arr(center_index) + arr(right_index)); /* Right */ + double bottom_value = -coeff3 * (att(center_index) + att(bottom_index)); /* Bottom */ + double top_value = -coeff4 * (att(center_index) + att(top_index)); /* Top */ + + double center_value = + (+0.25 * (h1 + h2) * (k1 + k2) * coeff_beta[center_index] * fabs(detDF(center_index)) /* beta_{i,j} */ + - left_value /* Center: (Left) */ + - right_value /* Center: (Right) */ + - bottom_value /* Center: (Bottom) */ + - top_value /* Center: (Top) */ + ); + + double bottom_left_value = -0.25 * (art(left_index) + art(bottom_index)); /* Bottom Left */ + double bottom_right_value = +0.25 * (art(right_index) + art(bottom_index)); /* Bottom Right */ + double top_left_value = +0.25 * (art(left_index) + art(top_index)); /* Top Left */ + double top_right_value = -0.25 * (art(right_index) + art(top_index)); /* Top Right */ + + /* 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 = center_value; + 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 = left_value; + updateMatrixElement(solver_matrix, ptr, offset, row, col, val); + } + + offset = CenterStencil[StencilPosition::Right]; + col = right_index; + val = right_value; + updateMatrixElement(solver_matrix, ptr, offset, row, col, val); + + offset = CenterStencil[StencilPosition::Bottom]; + col = bottom_index; + val = bottom_value; + updateMatrixElement(solver_matrix, ptr, offset, row, col, val); + + offset = CenterStencil[StencilPosition::Top]; + col = top_index; + val = top_value; + 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::BottomLeft]; + col = bottom_left_index; + val = bottom_left_value; + updateMatrixElement(solver_matrix, ptr, offset, row, col, val); + } + + offset = CenterStencil[StencilPosition::BottomRight]; + col = bottom_right_index; + val = bottom_right_value; + 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::TopLeft]; + col = top_left_index; + val = top_left_value; + updateMatrixElement(solver_matrix, ptr, offset, row, col, val); + } + + offset = CenterStencil[StencilPosition::TopRight]; + col = top_right_index; + val = top_right_value; + updateMatrixElement(solver_matrix, ptr, offset, row, col, val); + } + /* ------------------------------- */ + /* Node next to the outer boundary */ + /* ------------------------------- */ + else if (i_r == grid.nr() - 2) { + int i_theta_M1 = grid.wrapThetaIndex(i_theta - 1); + int i_theta_P1 = grid.wrapThetaIndex(i_theta + 1); + + double h1 = grid.radialSpacing(i_r - 1); + 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; + + int center_nz_index = getSolverMatrixIndex(i_r, i_theta); + + int center_index = grid.index(i_r, i_theta); + int left_index = grid.index(i_r - 1, i_theta); + int right_index = grid.index(i_r + 1, i_theta); + int bottom_index = grid.index(i_r, i_theta_M1); + int top_index = grid.index(i_r, i_theta_P1); + int bottom_left_index = grid.index(i_r - 1, i_theta_M1); + int bottom_right_index = grid.index(i_r + 1, i_theta_M1); + int top_left_index = grid.index(i_r - 1, i_theta_P1); + int top_right_index = grid.index(i_r + 1, i_theta_P1); + + double left_value = -coeff1 * (arr(center_index) + arr(left_index)); /* Left */ + double right_value = -coeff2 * (arr(center_index) + arr(right_index)); /* Right */ + double bottom_value = -coeff3 * (att(center_index) + att(bottom_index)); /* Bottom */ + double top_value = -coeff4 * (att(center_index) + att(top_index)); /* Top */ + + double center_value = + (+0.25 * (h1 + h2) * (k1 + k2) * coeff_beta[center_index] * fabs(detDF(center_index)) /* beta_{i,j} */ + - left_value /* Center: (Left) */ + - right_value /* Center: (Right) */ + - bottom_value /* Center: (Bottom) */ + - top_value /* Center: (Top) */ + ); + + double bottom_left_value = -0.25 * (art(left_index) + art(bottom_index)); /* Bottom Left */ + double bottom_right_value = +0.25 * (art(right_index) + art(bottom_index)); /* Bottom Right */ + double top_left_value = +0.25 * (art(left_index) + art(top_index)); /* Top Left */ + double top_right_value = -0.25 * (art(right_index) + art(top_index)); /* Top Right */ + + /* 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 = center_value; + updateMatrixElement(solver_matrix, ptr, offset, row, col, val); + + offset = CenterStencil[StencilPosition::Left]; + col = left_index; + val = left_value; + 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 = bottom_value; + updateMatrixElement(solver_matrix, ptr, offset, row, col, val); + + offset = CenterStencil[StencilPosition::Top]; + col = top_index; + val = top_value; + updateMatrixElement(solver_matrix, ptr, offset, row, col, val); + + offset = CenterStencil[StencilPosition::BottomLeft]; + col = bottom_left_index; + val = bottom_left_value; + updateMatrixElement(solver_matrix, ptr, offset, row, col, val); + + /* BottomRight REMOVED: Moved to the right hand side to make the matrix symmetric */ + + offset = CenterStencil[StencilPosition::TopLeft]; + col = top_left_index; + val = top_left_value; + updateMatrixElement(solver_matrix, ptr, offset, row, col, val); + + /* TopRight REMOVED: Moved to the right hand side to make the matrix symmetric */ + } + /* ------------------------------------ */ + /* Node on the outer dirichlet boundary */ + /* ------------------------------------ */ + else if (i_r == grid.nr() - 1) { + int center_nz_index = getSolverMatrixIndex(i_r, i_theta); + + int center_index = grid.index(i_r, i_theta); + + /* 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); + } +} void DirectSolver_COO_MUMPS_Take::buildSolverMatrixCircleSection(const int i_r, SparseMatrixCOO& solver_matrix) { assert(level_cache_.cacheDensityProfileCoefficients()); assert(level_cache_.cacheDomainGeometry()); - const auto& arr = level_cache_.arr(); - const auto& att = level_cache_.att(); - const auto& art = level_cache_.art(); - const auto& detDF = level_cache_.detDF(); - const auto& coeff_beta = level_cache_.coeff_beta(); + ConstVector arr = level_cache_.arr(); + ConstVector att = level_cache_.att(); + ConstVector art = level_cache_.art(); + ConstVector detDF = level_cache_.detDF(); + ConstVector coeff_beta = level_cache_.coeff_beta(); for (int i_theta = 0; i_theta < grid_.ntheta(); i_theta++) { // Build solver matrix at the current node - NODE_BUILD_SOLVER_MATRIX_TAKE(i_r, i_theta, grid_, DirBC_Interior_, solver_matrix, arr, att, art, detDF, - coeff_beta); + nodeBuildSolverMatrixTake(i_r, i_theta, grid_, DirBC_Interior_, solver_matrix, arr, att, art, detDF, + coeff_beta); } } @@ -471,16 +475,16 @@ void DirectSolver_COO_MUMPS_Take::buildSolverMatrixRadialSection(const int i_the assert(level_cache_.cacheDensityProfileCoefficients()); assert(level_cache_.cacheDomainGeometry()); - const auto& arr = level_cache_.arr(); - const auto& att = level_cache_.att(); - const auto& art = level_cache_.art(); - const auto& detDF = level_cache_.detDF(); - const auto& coeff_beta = level_cache_.coeff_beta(); + ConstVector arr = level_cache_.arr(); + ConstVector att = level_cache_.att(); + ConstVector art = level_cache_.art(); + ConstVector detDF = level_cache_.detDF(); + ConstVector coeff_beta = level_cache_.coeff_beta(); for (int i_r = grid_.numberSmootherCircles(); i_r < grid_.nr(); i_r++) { // Build solver matrix at the current node - NODE_BUILD_SOLVER_MATRIX_TAKE(i_r, i_theta, grid_, DirBC_Interior_, solver_matrix, arr, att, art, detDF, - coeff_beta); + nodeBuildSolverMatrixTake(i_r, i_theta, grid_, DirBC_Interior_, solver_matrix, arr, att, art, detDF, + coeff_beta); } } diff --git a/src/DirectSolver/DirectSolver-CSR-LU-Take/buildSolverMatrix.cpp b/src/DirectSolver/DirectSolver-CSR-LU-Take/buildSolverMatrix.cpp index efef52c3..efd03176 100644 --- a/src/DirectSolver/DirectSolver-CSR-LU-Take/buildSolverMatrix.cpp +++ b/src/DirectSolver/DirectSolver-CSR-LU-Take/buildSolverMatrix.cpp @@ -1,256 +1,260 @@ #include "../../../include/DirectSolver/DirectSolver-CSR-LU-Take/directSolverTakeCustomLU.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_TAKE(i_r, i_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 > 0 && i_r < grid.nr() - 1) { \ - 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); \ - const int bottom_left_index = grid.index(i_r - 1, i_theta_M1); \ - const int bottom_right_index = grid.index(i_r + 1, i_theta_M1); \ - const int top_left_index = grid.index(i_r - 1, i_theta_P1); \ - const int top_right_index = grid.index(i_r + 1, i_theta_P1); \ - \ - const double left_value = -coeff1 * (arr[center_index] + arr[left_index]); /* Left */ \ - const double right_value = -coeff2 * (arr[center_index] + arr[right_index]); /* Right */ \ - const double bottom_value = -coeff3 * (att[center_index] + att[bottom_index]); /* Bottom */ \ - const double top_value = -coeff4 * (att[center_index] + att[top_index]); /* Top */ \ - \ - const double center_value = \ - (+0.25 * (h1 + h2) * (k1 + k2) * coeff_beta[center_index] * fabs(detDF[center_index]) /* beta_{i,j} */ \ - - left_value /* Center: (Left) */ \ - - right_value /* Center: (Right) */ \ - - bottom_value /* Center: (Bottom) */ \ - - top_value /* Center: (Top) */ \ - ); \ - \ - const double bottom_left_value = -0.25 * (art[left_index] + art[bottom_index]); /* Bottom Left */ \ - const double bottom_right_value = +0.25 * (art[right_index] + art[bottom_index]); /* Bottom Right */ \ - const double top_left_value = +0.25 * (art[left_index] + art[top_index]); /* Top Left */ \ - const double top_right_value = -0.25 * (art[right_index] + art[top_index]); /* Top Right */ \ - \ - /* Fill matrix row of (i,j) */ \ - row = center_index; \ - \ - const Stencil& CenterStencil = getStencil(i_r); \ - \ - offset = CenterStencil[StencilPosition::Center]; \ - col = center_index; \ - val = center_value; \ - UPDATE_MATRIX_ELEMENT(solver_matrix, offset, row, col, val); \ - \ - offset = CenterStencil[StencilPosition::Left]; \ - col = left_index; \ - val = left_value; \ - UPDATE_MATRIX_ELEMENT(solver_matrix, offset, row, col, val); \ - \ - offset = CenterStencil[StencilPosition::Right]; \ - col = right_index; \ - val = right_value; \ - UPDATE_MATRIX_ELEMENT(solver_matrix, offset, row, col, val); \ - \ - offset = CenterStencil[StencilPosition::Bottom]; \ - col = bottom_index; \ - val = bottom_value; \ - UPDATE_MATRIX_ELEMENT(solver_matrix, offset, row, col, val); \ - \ - offset = CenterStencil[StencilPosition::Top]; \ - col = top_index; \ - val = top_value; \ - UPDATE_MATRIX_ELEMENT(solver_matrix, offset, row, col, val); \ - \ - offset = CenterStencil[StencilPosition::BottomLeft]; \ - col = bottom_left_index; \ - val = bottom_left_value; \ - UPDATE_MATRIX_ELEMENT(solver_matrix, offset, row, col, val); \ - \ - offset = CenterStencil[StencilPosition::BottomRight]; \ - col = bottom_right_index; \ - val = bottom_right_value; \ - UPDATE_MATRIX_ELEMENT(solver_matrix, offset, row, col, val); \ - \ - offset = CenterStencil[StencilPosition::TopLeft]; \ - col = top_left_index; \ - val = top_left_value; \ - UPDATE_MATRIX_ELEMENT(solver_matrix, offset, row, col, val); \ - \ - offset = CenterStencil[StencilPosition::TopRight]; \ - col = top_right_index; \ - val = top_right_value; \ - 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 int center_index = grid.index(i_r, i_theta); \ - \ - /* 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); \ - } \ - 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); \ - \ - 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, 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); \ - const int bottom_right_index = grid.index(i_r + 1, i_theta_M1); \ - const int top_right_index = grid.index(i_r + 1, i_theta_P1); \ - \ - const double left_value = -coeff1 * (arr[center_index] + arr[left_index]); /* Left */ \ - const double right_value = -coeff2 * (arr[center_index] + arr[right_index]); /* Right */ \ - const double bottom_value = -coeff3 * (att[center_index] + att[bottom_index]); /* Bottom */ \ - const double top_value = -coeff4 * (att[center_index] + att[top_index]); /* Top */ \ - \ - const double center_value = (+0.25 * (h1 + h2) * (k1 + k2) * coeff_beta[center_index] * \ - fabs(detDF[center_index]) /* beta_{i,j} */ \ - - left_value /* Center: (Left) */ \ - - right_value /* Center: (Right) */ \ - - bottom_value /* Center: (Bottom) */ \ - - top_value /* Center: (Top) */ \ - ); \ - \ - const double bottom_right_value = +0.25 * (art[right_index] + art[bottom_index]); /* Bottom Right */ \ - const double top_right_value = -0.25 * (art[right_index] + art[top_index]); /* Top Right */ \ - \ - /* Fill matrix row of (i,j) */ \ - row = center_index; \ - \ - const Stencil& CenterStencil = getStencil(i_r); \ - \ - offset = CenterStencil[StencilPosition::Center]; \ - col = center_index; \ - val = center_value; \ - UPDATE_MATRIX_ELEMENT(solver_matrix, offset, row, col, val); \ - \ - offset = CenterStencil[StencilPosition::Left]; \ - col = left_index; \ - val = left_value; \ - UPDATE_MATRIX_ELEMENT(solver_matrix, offset, row, col, val); \ - \ - offset = CenterStencil[StencilPosition::Right]; \ - col = right_index; \ - val = right_value; \ - UPDATE_MATRIX_ELEMENT(solver_matrix, offset, row, col, val); \ - \ - offset = CenterStencil[StencilPosition::Bottom]; \ - col = bottom_index; \ - val = bottom_value; \ - UPDATE_MATRIX_ELEMENT(solver_matrix, offset, row, col, val); \ - \ - offset = CenterStencil[StencilPosition::Top]; \ - col = top_index; \ - val = top_value; \ - UPDATE_MATRIX_ELEMENT(solver_matrix, offset, row, col, val); \ - \ - /* BottomLeft: REMOVED DUE TO ARTIFICAL 7 POINT STENCIL */ \ - \ - offset = CenterStencil[StencilPosition::BottomRight]; \ - col = bottom_right_index; \ - val = bottom_right_value; \ - UPDATE_MATRIX_ELEMENT(solver_matrix, offset, row, col, val); \ - \ - /* TopLeft: REMOVED DUE TO ARTIFICAL 7 POINT STENCIL */ \ - \ - offset = CenterStencil[StencilPosition::TopRight]; \ - col = top_right_index; \ - val = top_right_value; \ - UPDATE_MATRIX_ELEMENT(solver_matrix, offset, row, col, val); \ - } \ - } \ - /* ------------------------------------ */ \ - /* Node on the outer dirichlet boundary */ \ - /* ------------------------------------ */ \ - else if (i_r == grid.nr() - 1) { \ - const int center_index = grid.index(i_r, i_theta); \ - \ - /* 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); \ - } \ - } 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_Take::nodeBuildSolverMatrixTake(int i_r, int i_theta, const PolarGrid& grid, + bool DirBC_Interior, SparseMatrixCSR& solver_matrix, + ConstVector& arr, ConstVector& att, + ConstVector& art, ConstVector& detDF, + ConstVector& coeff_beta) +{ + int offset; + int row, col; + double val; + /* -------------------- */ + /* Node in the interior */ + /* -------------------- */ + if (i_r > 0 && i_r < grid.nr() - 1) { + int i_theta_M1 = grid.wrapThetaIndex(i_theta - 1); + int i_theta_P1 = grid.wrapThetaIndex(i_theta + 1); + + double h1 = grid.radialSpacing(i_r - 1); + 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; + + int center_index = grid.index(i_r, i_theta); + int left_index = grid.index(i_r - 1, i_theta); + int right_index = grid.index(i_r + 1, i_theta); + int bottom_index = grid.index(i_r, i_theta_M1); + int top_index = grid.index(i_r, i_theta_P1); + int bottom_left_index = grid.index(i_r - 1, i_theta_M1); + int bottom_right_index = grid.index(i_r + 1, i_theta_M1); + int top_left_index = grid.index(i_r - 1, i_theta_P1); + int top_right_index = grid.index(i_r + 1, i_theta_P1); + + double left_value = -coeff1 * (arr[center_index] + arr[left_index]); /* Left */ + double right_value = -coeff2 * (arr[center_index] + arr[right_index]); /* Right */ + double bottom_value = -coeff3 * (att[center_index] + att[bottom_index]); /* Bottom */ + double top_value = -coeff4 * (att[center_index] + att[top_index]); /* Top */ + + double center_value = + (+0.25 * (h1 + h2) * (k1 + k2) * coeff_beta[center_index] * std::fabs(detDF[center_index]) /* beta_{i,j} */ + - left_value /* Center: (Left) */ + - right_value /* Center: (Right) */ + - bottom_value /* Center: (Bottom) */ + - top_value /* Center: (Top) */ + ); + + double bottom_left_value = -0.25 * (art[left_index] + art[bottom_index]); /* Bottom Left */ + double bottom_right_value = +0.25 * (art[right_index] + art[bottom_index]); /* Bottom Right */ + double top_left_value = +0.25 * (art[left_index] + art[top_index]); /* Top Left */ + double top_right_value = -0.25 * (art[right_index] + art[top_index]); /* Top Right */ + + /* Fill matrix row of (i,j) */ + row = center_index; + + const Stencil& CenterStencil = getStencil(i_r); + + offset = CenterStencil[StencilPosition::Center]; + col = center_index; + val = center_value; + updateMatrixElement(solver_matrix, offset, row, col, val); + + offset = CenterStencil[StencilPosition::Left]; + col = left_index; + val = left_value; + updateMatrixElement(solver_matrix, offset, row, col, val); + + offset = CenterStencil[StencilPosition::Right]; + col = right_index; + val = right_value; + updateMatrixElement(solver_matrix, offset, row, col, val); + + offset = CenterStencil[StencilPosition::Bottom]; + col = bottom_index; + val = bottom_value; + updateMatrixElement(solver_matrix, offset, row, col, val); + + offset = CenterStencil[StencilPosition::Top]; + col = top_index; + val = top_value; + updateMatrixElement(solver_matrix, offset, row, col, val); + + offset = CenterStencil[StencilPosition::BottomLeft]; + col = bottom_left_index; + val = bottom_left_value; + updateMatrixElement(solver_matrix, offset, row, col, val); + + offset = CenterStencil[StencilPosition::BottomRight]; + col = bottom_right_index; + val = bottom_right_value; + updateMatrixElement(solver_matrix, offset, row, col, val); + + offset = CenterStencil[StencilPosition::TopLeft]; + col = top_left_index; + val = top_left_value; + updateMatrixElement(solver_matrix, offset, row, col, val); + + offset = CenterStencil[StencilPosition::TopRight]; + col = top_right_index; + val = top_right_value; + 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) { + int center_index = grid.index(i_r, i_theta); + + /* 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); + } + 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. + assert(grid_.ntheta() % 2 == 0); + + int i_theta_M1 = grid.wrapThetaIndex(i_theta - 1); + int i_theta_P1 = grid.wrapThetaIndex(i_theta + 1); + 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; + + int center_index = grid.index(i_r, i_theta); + int left_index = grid.index(i_r, i_theta_AcrossOrigin); + int right_index = grid.index(i_r + 1, i_theta); + int bottom_index = grid.index(i_r, i_theta_M1); + int top_index = grid.index(i_r, i_theta_P1); + int bottom_right_index = grid.index(i_r + 1, i_theta_M1); + int top_right_index = grid.index(i_r + 1, i_theta_P1); + + double left_value = -coeff1 * (arr[center_index] + arr[left_index]); /* Left */ + double right_value = -coeff2 * (arr[center_index] + arr[right_index]); /* Right */ + double bottom_value = -coeff3 * (att[center_index] + att[bottom_index]); /* Bottom */ + double top_value = -coeff4 * (att[center_index] + att[top_index]); /* Top */ + + double center_value = + (+0.25 * (h1 + h2) * (k1 + k2) * coeff_beta[center_index] * fabs(detDF[center_index]) /* beta_{i,j} */ + - left_value /* Center: (Left) */ + - right_value /* Center: (Right) */ + - bottom_value /* Center: (Bottom) */ + - top_value /* Center: (Top) */ + ); + + double bottom_right_value = +0.25 * (art[right_index] + art[bottom_index]); /* Bottom Right */ + double top_right_value = -0.25 * (art[right_index] + art[top_index]); /* Top Right */ + + /* Fill matrix row of (i,j) */ + row = center_index; + + const Stencil& CenterStencil = getStencil(i_r); + + offset = CenterStencil[StencilPosition::Center]; + col = center_index; + val = center_value; + updateMatrixElement(solver_matrix, offset, row, col, val); + + offset = CenterStencil[StencilPosition::Left]; + col = left_index; + val = left_value; + updateMatrixElement(solver_matrix, offset, row, col, val); + + offset = CenterStencil[StencilPosition::Right]; + col = right_index; + val = right_value; + updateMatrixElement(solver_matrix, offset, row, col, val); + + offset = CenterStencil[StencilPosition::Bottom]; + col = bottom_index; + val = bottom_value; + updateMatrixElement(solver_matrix, offset, row, col, val); + + offset = CenterStencil[StencilPosition::Top]; + col = top_index; + val = top_value; + updateMatrixElement(solver_matrix, offset, row, col, val); + + /* BottomLeft: REMOVED DUE TO ARTIFICAL 7 POINT STENCIL */ + + offset = CenterStencil[StencilPosition::BottomRight]; + col = bottom_right_index; + val = bottom_right_value; + updateMatrixElement(solver_matrix, offset, row, col, val); + + /* TopLeft: REMOVED DUE TO ARTIFICAL 7 POINT STENCIL */ + + offset = CenterStencil[StencilPosition::TopRight]; + col = top_right_index; + val = top_right_value; + updateMatrixElement(solver_matrix, offset, row, col, val); + } + } + /* ------------------------------------ */ + /* Node on the outer dirichlet boundary */ + /* ------------------------------------ */ + else if (i_r == grid.nr() - 1) { + int center_index = grid.index(i_r, i_theta); + + /* 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); + } +} void DirectSolver_CSR_LU_Take::buildSolverMatrixCircleSection(const int i_r, SparseMatrixCSR& solver_matrix) { assert(level_cache_.cacheDensityProfileCoefficients()); assert(level_cache_.cacheDomainGeometry()); - const auto& arr = level_cache_.arr(); - const auto& att = level_cache_.att(); - const auto& art = level_cache_.art(); - const auto& detDF = level_cache_.detDF(); - const auto& coeff_beta = level_cache_.coeff_beta(); + ConstVector arr = level_cache_.arr(); + ConstVector att = level_cache_.att(); + ConstVector art = level_cache_.art(); + ConstVector detDF = level_cache_.detDF(); + ConstVector coeff_beta = level_cache_.coeff_beta(); for (int i_theta = 0; i_theta < grid_.ntheta(); i_theta++) { // Build solver matrix at the current node - NODE_BUILD_SOLVER_MATRIX_TAKE(i_r, i_theta, grid_, DirBC_Interior_, solver_matrix, arr, att, art, detDF, - coeff_beta); + nodeBuildSolverMatrixTake(i_r, i_theta, grid_, DirBC_Interior_, solver_matrix, arr, att, art, detDF, + coeff_beta); } } @@ -259,16 +263,16 @@ void DirectSolver_CSR_LU_Take::buildSolverMatrixRadialSection(const int i_theta, assert(level_cache_.cacheDensityProfileCoefficients()); assert(level_cache_.cacheDomainGeometry()); - const auto& arr = level_cache_.arr(); - const auto& att = level_cache_.att(); - const auto& art = level_cache_.art(); - const auto& detDF = level_cache_.detDF(); - const auto& coeff_beta = level_cache_.coeff_beta(); + ConstVector arr = level_cache_.arr(); + ConstVector att = level_cache_.att(); + ConstVector art = level_cache_.art(); + ConstVector detDF = level_cache_.detDF(); + ConstVector coeff_beta = level_cache_.coeff_beta(); for (int i_r = grid_.numberSmootherCircles(); i_r < grid_.nr(); i_r++) { // Build solver matrix at the current node - NODE_BUILD_SOLVER_MATRIX_TAKE(i_r, i_theta, grid_, DirBC_Interior_, solver_matrix, arr, att, art, detDF, - coeff_beta); + nodeBuildSolverMatrixTake(i_r, i_theta, grid_, DirBC_Interior_, solver_matrix, arr, att, art, detDF, + coeff_beta); } }