diff --git a/include/DirectSolver/DirectSolverGive/buildSolverMatrix.inl b/include/DirectSolver/DirectSolverGive/buildSolverMatrix.inl index 59637462..48b2877a 100644 --- a/include/DirectSolver/DirectSolverGive/buildSolverMatrix.inl +++ b/include/DirectSolver/DirectSolverGive/buildSolverMatrix.inl @@ -5,8 +5,8 @@ namespace direct_solver_coo_mumps_give #ifdef GMGPOLAR_USE_MUMPS // When using the MUMPS solver, the matrix is assembled in COO format. -static inline void updateMatrixElement(SparseMatrixCOO& matrix, int ptr, int offset, int row, int column, - double value) +static inline void updateMatrixElement(SparseMatrixCOO& matrix, int ptr, int offset, int row, + int column, double value) { matrix.set_row_index(ptr + offset, row); matrix.set_col_index(ptr + offset, column); diff --git a/include/DirectSolver/DirectSolverTake/buildSolverMatrix.inl b/include/DirectSolver/DirectSolverTake/buildSolverMatrix.inl index 4cdb4487..a144a9c6 100644 --- a/include/DirectSolver/DirectSolverTake/buildSolverMatrix.inl +++ b/include/DirectSolver/DirectSolverTake/buildSolverMatrix.inl @@ -5,8 +5,8 @@ namespace direct_solver_coo_mumps_take #ifdef GMGPOLAR_USE_MUMPS // When using the MUMPS solver, the matrix is assembled in COO format. -static inline void updateMatrixElement(SparseMatrixCOO& matrix, int ptr, int offset, int row, int column, - double value) +static inline void updateMatrixElement(SparseMatrixCOO& matrix, int ptr, int offset, int row, + int column, double value) { matrix.set_row_index(ptr + offset, row); matrix.set_col_index(ptr + offset, column); diff --git a/include/ExtrapolatedSmoother/ExtrapolatedSmootherGive/buildInnerBoundaryAsc.inl b/include/ExtrapolatedSmoother/ExtrapolatedSmootherGive/buildInnerBoundaryAsc.inl index 0ff8b1b2..ad6b3362 100644 --- a/include/ExtrapolatedSmoother/ExtrapolatedSmootherGive/buildInnerBoundaryAsc.inl +++ b/include/ExtrapolatedSmoother/ExtrapolatedSmootherGive/buildInnerBoundaryAsc.inl @@ -7,8 +7,8 @@ namespace extrapolated_smoother_give #ifdef GMGPOLAR_USE_MUMPS // When using the MUMPS solver, the matrix is assembled in COO format. -static inline void update_CSR_COO_MatrixElement(SparseMatrixCOO& matrix, int ptr, int offset, int row, - int column, double value) +static inline void update_CSR_COO_MatrixElement(SparseMatrixCOO& matrix, int ptr, int offset, + int row, int column, double value) { matrix.set_row_index(ptr + offset, row); matrix.set_col_index(ptr + offset, column); diff --git a/include/ExtrapolatedSmoother/ExtrapolatedSmootherTake/buildInnerBoundaryAsc.inl b/include/ExtrapolatedSmoother/ExtrapolatedSmootherTake/buildInnerBoundaryAsc.inl index aae19b48..ebf2eb4d 100644 --- a/include/ExtrapolatedSmoother/ExtrapolatedSmootherTake/buildInnerBoundaryAsc.inl +++ b/include/ExtrapolatedSmoother/ExtrapolatedSmootherTake/buildInnerBoundaryAsc.inl @@ -5,8 +5,8 @@ namespace extrapolated_smoother_take #ifdef GMGPOLAR_USE_MUMPS // When using the MUMPS solver, the matrix is assembled in COO format. -static inline void update_CSR_COO_MatrixElement(SparseMatrixCOO& matrix, int ptr, int offset, int row, - int column, double value) +static inline void update_CSR_COO_MatrixElement(SparseMatrixCOO& matrix, int ptr, int offset, + int row, int column, double value) { matrix.set_row_index(ptr + offset, row); matrix.set_col_index(ptr + offset, column); diff --git a/include/LinearAlgebra/Solvers/tridiagonal_solver.h b/include/LinearAlgebra/Solvers/tridiagonal_solver.h index 77d8601c..ce98eb8c 100644 --- a/include/LinearAlgebra/Solvers/tridiagonal_solver.h +++ b/include/LinearAlgebra/Solvers/tridiagonal_solver.h @@ -54,6 +54,16 @@ class BatchedTridiagonalSolver { return main_diagonal_(batch_idx * matrix_dimension_ + index); } + KOKKOS_INLINE_FUNCTION + void set_main_diagonal(const int batch_idx, const int index, const T value) const + { + main_diagonal_(batch_idx * matrix_dimension_ + index) = value; + } + KOKKOS_INLINE_FUNCTION + void add_main_diagonal(const int batch_idx, const int index, const T value) const + { + main_diagonal_(batch_idx * matrix_dimension_ + index) += value; + } KOKKOS_INLINE_FUNCTION const T& sub_diagonal(const int batch_idx, const int index) const @@ -65,18 +75,37 @@ class BatchedTridiagonalSolver { return sub_diagonal_(batch_idx * matrix_dimension_ + index); } + KOKKOS_INLINE_FUNCTION + void set_sub_diagonal(const int batch_idx, const int index, const T value) const + { + sub_diagonal_(batch_idx * matrix_dimension_ + index) = value; + } + KOKKOS_INLINE_FUNCTION + void add_sub_diagonal(const int batch_idx, const int index, const T value) const + { + sub_diagonal_(batch_idx * matrix_dimension_ + index) += value; + } KOKKOS_INLINE_FUNCTION const T& cyclic_corner(const int batch_idx) const { return sub_diagonal_(batch_idx * matrix_dimension_ + (matrix_dimension_ - 1)); } - KOKKOS_INLINE_FUNCTION T& cyclic_corner(const int batch_idx) { return sub_diagonal_(batch_idx * matrix_dimension_ + (matrix_dimension_ - 1)); } + KOKKOS_INLINE_FUNCTION + void set_cyclic_corner(const int batch_idx, const T value) const + { + sub_diagonal_(batch_idx * matrix_dimension_ + (matrix_dimension_ - 1)) = value; + } + KOKKOS_INLINE_FUNCTION + void add_cyclic_corner(const int batch_idx, const T value) const + { + sub_diagonal_(batch_idx * matrix_dimension_ + (matrix_dimension_ - 1)) += value; + } /* ---------------------------------------------- */ /* Setup: Cholesky Decomposition: A = L * D * L^T */ @@ -233,6 +262,7 @@ class BatchedTridiagonalSolver rhs(offset + 0) + cyclic_corner_element / gamma(batch_idx) * rhs(offset + matrix_dimension - 1); const T dot_product_u_v = buffer(offset + 0) + cyclic_corner_element / gamma(batch_idx) * buffer(offset + matrix_dimension - 1); + const T factor = dot_product_x_v / (1.0 + dot_product_u_v); for (int i = 0; i < matrix_dimension; i++) { @@ -305,10 +335,10 @@ class BatchedTridiagonalSolver int matrix_dimension_; int batch_count_; - Vector main_diagonal_; - Vector sub_diagonal_; - Vector buffer_; - Vector gamma_; + AllocatableVector main_diagonal_; + AllocatableVector sub_diagonal_; + AllocatableVector buffer_; + AllocatableVector gamma_; bool is_cyclic_; bool is_factorized_; diff --git a/include/Smoother/SmootherGive/buildInnerBoundaryAsc.inl b/include/Smoother/SmootherGive/buildInnerBoundaryAsc.inl index dd92f089..c59c0359 100644 --- a/include/Smoother/SmootherGive/buildInnerBoundaryAsc.inl +++ b/include/Smoother/SmootherGive/buildInnerBoundaryAsc.inl @@ -5,8 +5,8 @@ namespace smoother_give #ifdef GMGPOLAR_USE_MUMPS // When using the MUMPS solver, the matrix is assembled in COO format. -static inline void update_CSR_COO_MatrixElement(SparseMatrixCOO& matrix, int ptr, int offset, int row, - int column, double value) +static inline void update_CSR_COO_MatrixElement(SparseMatrixCOO& matrix, int ptr, int offset, + int row, int column, double value) { matrix.set_row_index(ptr + offset, row); matrix.set_col_index(ptr + offset, column); diff --git a/include/Smoother/SmootherTake/buildInnerBoundaryAsc.inl b/include/Smoother/SmootherTake/buildInnerBoundaryAsc.inl index b2f3cfd3..60fcf88b 100644 --- a/include/Smoother/SmootherTake/buildInnerBoundaryAsc.inl +++ b/include/Smoother/SmootherTake/buildInnerBoundaryAsc.inl @@ -2,33 +2,35 @@ namespace smoother_take { +#include "matrixStencil.inl" #ifdef GMGPOLAR_USE_MUMPS // When using the MUMPS solver, the matrix is assembled in COO format. -static inline void update_CSR_COO_MatrixElement(SparseMatrixCOO& matrix, int ptr, int offset, int row, - int column, double value) +static inline void update_CSR_COO_MatrixElement(const SparseMatrixCOO& matrix, int ptr, + int offset, int row, int column, double value) { matrix.set_row_index(ptr + offset, row); matrix.set_col_index(ptr + offset, column); - matrix.set_value(ptr + offset , value); + matrix.set_value(ptr + offset, value); } #else // When using the in-house solver, the matrix is stored in CSR format. -static inline void update_CSR_COO_MatrixElement(SparseMatrixCSR& matrix, int ptr, int offset, - int row, int column, double value) +static inline void update_CSR_COO_MatrixElement(const SparseMatrixCSR& matrix, int ptr, + int offset, int row, int column, double value) { matrix.set_row_nz_index(row, offset, column); matrix.set_row_nz_entry(row, offset, value); } #endif -} // namespace smoother_take - -template -void SmootherTake::nodeBuildInteriorBoundarySolverMatrix( - int i_theta, const PolarGrid& grid, bool DirBC_Interior, InnerBoundaryMatrix& matrix, ConstVector& arr, - ConstVector& att, ConstVector& art, ConstVector& detDF, ConstVector& coeff_beta) +template +static inline void nodeBuildInteriorBoundarySolverMatrix(int i_theta, const PolarGrid& grid, bool DirBC_Interior, + const MatrixType& matrix, ConstVector& arr, + ConstVector& att, ConstVector& art, + ConstVector& detDF, ConstVector& coeff_beta) { + using smoother_take::getCircleAscIndex; + using smoother_take::getStencil; using smoother_take::update_CSR_COO_MatrixElement; assert(i_theta >= 0 && i_theta < grid.ntheta()); @@ -47,13 +49,13 @@ void SmootherTake::nodeBuildInteriorBoundarySolverMatrix( /* ------------------------------------------------ */ if (DirBC_Interior) { const int center_index = i_theta; - const int center_nz_index = getCircleAscIndex(i_r, i_theta); + const int center_nz_index = getCircleAscIndex(i_r, i_theta, DirBC_Interior); /* Fill matrix row of (i,j) */ row = center_index; ptr = center_nz_index; - const Stencil& CenterStencil = getStencil(i_r); + const Stencil& CenterStencil = getStencil(i_r, DirBC_Interior); offset = CenterStencil[StencilPosition::Center]; column = center_index; @@ -92,7 +94,7 @@ void SmootherTake::nodeBuildInteriorBoundarySolverMatrix( const int bottom_index = i_theta_M1; const int top_index = i_theta_P1; - const int center_nz_index = getCircleAscIndex(i_r, i_theta); + const int center_nz_index = getCircleAscIndex(i_r, i_theta, DirBC_Interior); const double left_value = -coeff1 * (arr[center] + arr[left]); const double right_value = -coeff2 * (arr[center] + arr[right]); @@ -106,7 +108,7 @@ void SmootherTake::nodeBuildInteriorBoundarySolverMatrix( row = center_index; ptr = center_nz_index; - const Stencil& CenterStencil = getStencil(i_r); + const Stencil& CenterStencil = getStencil(i_r, DirBC_Interior); offset = CenterStencil[StencilPosition::Center]; column = center_index; @@ -130,10 +132,14 @@ void SmootherTake::nodeBuildInteriorBoundarySolverMatrix( } } +} // namespace smoother_take + template typename SmootherTake::InnerBoundaryMatrix SmootherTake::buildInteriorBoundarySolverMatrix() { + using smoother_take::nodeBuildInteriorBoundarySolverMatrix; + const PolarGrid& grid = Smoother::grid_; const LevelCacheType& level_cache = Smoother::level_cache_; const bool DirBC_Interior = Smoother::DirBC_Interior_; @@ -149,7 +155,8 @@ SmootherTake::buildInteriorBoundarySolverMatrix() // the COO_Mumps_Solver optimizes the factorization by only using the upper triangular part of the matrix, // which is extracted by the COO_Mumps_Solver internally. #ifdef GMGPOLAR_USE_MUMPS - const int nnz = getNonZeroCountCircleAsc(i_r); + using smoother_take::getNonZeroCountCircleAsc; + const int nnz = getNonZeroCountCircleAsc(i_r, grid, DirBC_Interior); SparseMatrixCOO inner_boundary_solver_matrix(ntheta, ntheta, nnz); inner_boundary_solver_matrix.is_symmetric(true); #else @@ -169,11 +176,18 @@ SmootherTake::buildInteriorBoundarySolverMatrix() ConstVector detDF = level_cache.detDF(); ConstVector coeff_beta = level_cache.coeff_beta(); -#pragma omp parallel for num_threads(num_omp_threads) - for (int i_theta = 0; i_theta < ntheta; i_theta++) { - nodeBuildInteriorBoundarySolverMatrix(i_theta, grid, DirBC_Interior, inner_boundary_solver_matrix, arr, att, - art, detDF, coeff_beta); - } + const PolarGrid* grid_ptr = &grid; + + // Capture pointer to matrix for use in Kokkos parallel loop. + const InnerBoundaryMatrix* matrix_ptr = &inner_boundary_solver_matrix; + + Kokkos::parallel_for( + "Smoother Take: BuildInnerBoundaryAsc", Kokkos::RangePolicy<>(0, ntheta), KOKKOS_LAMBDA(const int i_theta) { + nodeBuildInteriorBoundarySolverMatrix(i_theta, *grid_ptr, DirBC_Interior, *matrix_ptr, arr, att, art, detDF, + coeff_beta); + }); + + Kokkos::fence(); return inner_boundary_solver_matrix; } diff --git a/include/Smoother/SmootherTake/buildTridiagonalAsc.inl b/include/Smoother/SmootherTake/buildTridiagonalAsc.inl index 9fd360b2..1867e79c 100644 --- a/include/Smoother/SmootherTake/buildTridiagonalAsc.inl +++ b/include/Smoother/SmootherTake/buildTridiagonalAsc.inl @@ -3,25 +3,24 @@ namespace smoother_take { -static inline void updateMatrixElement(BatchedTridiagonalSolver& solver, int batch, int row, int column, +static inline void updateMatrixElement(const BatchedTridiagonalSolver& solver, int batch, int row, int column, double value) { if (row == column) - solver.main_diagonal(batch, row) = value; + solver.set_main_diagonal(batch, row, value); else if (row == column - 1) - solver.sub_diagonal(batch, row) = value; + solver.set_sub_diagonal(batch, row, value); else if (row == 0 && column == solver.matrixDimension() - 1) - solver.cyclic_corner(batch) = value; + solver.set_cyclic_corner(batch, value); } -} // namespace smoother_take - -template -void SmootherTake::nodeBuildTridiagonalSolverMatrices( - int i_r, int i_theta, const PolarGrid& grid, bool DirBC_Interior, - BatchedTridiagonalSolver& circle_tridiagonal_solver, - BatchedTridiagonalSolver& radial_tridiagonal_solver, ConstVector& arr, ConstVector& att, - ConstVector& art, ConstVector& detDF, ConstVector& coeff_beta) +// Build the tridiagonal solver matrices for a specific node (i_r, i_theta) +static inline void nodeBuildTridiagonalSolverMatrices(int i_r, int i_theta, const PolarGrid& grid, bool DirBC_Interior, + const BatchedTridiagonalSolver& circle_tridiagonal_solver, + const BatchedTridiagonalSolver& radial_tridiagonal_solver, + ConstVector& arr, ConstVector& att, + ConstVector& art, ConstVector& detDF, + ConstVector& coeff_beta) { using smoother_take::updateMatrixElement; @@ -302,14 +301,21 @@ void SmootherTake::nodeBuildTridiagonalSolverMatrices( } } +} // namespace smoother_take + template void SmootherTake::buildTridiagonalSolverMatrices() { + using smoother_take::nodeBuildTridiagonalSolverMatrices; + const PolarGrid& grid = Smoother::grid_; const LevelCacheType& level_cache = Smoother::level_cache_; const bool DirBC_Interior = Smoother::DirBC_Interior_; const int num_omp_threads = Smoother::num_omp_threads_; + const BatchedTridiagonalSolver* circle_tridiagonal_solver_ptr = &circle_tridiagonal_solver_; + const BatchedTridiagonalSolver* radial_tridiagonal_solver_ptr = &radial_tridiagonal_solver_; + assert(level_cache.cacheDensityProfileCoefficients()); assert(level_cache.cacheDomainGeometry()); @@ -319,22 +325,36 @@ void SmootherTake::buildTridiagonalSolverMatrices() ConstVector detDF = level_cache.detDF(); ConstVector coeff_beta = level_cache.coeff_beta(); -#pragma omp parallel num_threads(num_omp_threads) - { -#pragma omp for nowait - for (int i_r = 0; i_r < grid.numberSmootherCircles(); i_r++) { - for (int i_theta = 0; i_theta < grid.ntheta(); i_theta++) { - nodeBuildTridiagonalSolverMatrices(i_r, i_theta, grid, DirBC_Interior, circle_tridiagonal_solver_, - radial_tridiagonal_solver_, arr, att, art, detDF, coeff_beta); - } - } - -#pragma omp for nowait - for (int i_theta = 0; i_theta < grid.ntheta(); i_theta++) { - for (int i_r = grid.numberSmootherCircles(); i_r < grid.nr(); i_r++) { - nodeBuildTridiagonalSolverMatrices(i_r, i_theta, grid, DirBC_Interior, circle_tridiagonal_solver_, - radial_tridiagonal_solver_, arr, att, art, detDF, coeff_beta); - } - } - } + const PolarGrid* grid_ptr = &grid; + + /* We split the loops into two regions to better respect the */ + /* access patterns of the smoother and improve cache locality. */ + + // The For loop matches circular access pattern */ + Kokkos::parallel_for( + "Smoother Take: BuildTridiagonalAsc (Circular)", + Kokkos::MDRangePolicy>( // Rank of the index space + {0, 0}, // Starting point of the index space + {grid.numberSmootherCircles(), grid.ntheta()} // Ending point of the index space + ), + // Kokkos lambda function to execute for each point in the index space + KOKKOS_LAMBDA(const int i_r, const int i_theta) { + nodeBuildTridiagonalSolverMatrices(i_r, i_theta, *grid_ptr, DirBC_Interior, circle_tridiagonal_solver_, + radial_tridiagonal_solver_, arr, att, art, detDF, coeff_beta); + }); + + /* For loop matches radial access pattern */ + Kokkos::parallel_for( + "Smoother Take: BuildTridiagonalAsc (Radial)", + Kokkos::MDRangePolicy>( // Rank of the index space + {0, grid.numberSmootherCircles()}, // Starting point of the index space + {grid.ntheta(), grid.nr()} // Ending point of the index space + ), + // Kokkos lambda function to execute for each point in the index space + KOKKOS_LAMBDA(const int i_theta, const int i_r) { + nodeBuildTridiagonalSolverMatrices(i_r, i_theta, *grid_ptr, DirBC_Interior, circle_tridiagonal_solver_, + radial_tridiagonal_solver_, arr, att, art, detDF, coeff_beta); + }); + + Kokkos::fence(); } \ No newline at end of file diff --git a/include/Smoother/SmootherTake/matrixStencil.inl b/include/Smoother/SmootherTake/matrixStencil.inl index a2135d66..bc11625f 100644 --- a/include/Smoother/SmootherTake/matrixStencil.inl +++ b/include/Smoother/SmootherTake/matrixStencil.inl @@ -1,7 +1,9 @@ -#pragma once +/* -------------- */ +/* Stencil access */ +/* -------------- */ -template -const Stencil& SmootherTake::getStencil(int i_r) const +// Select correct stencil depending on the grid position. +static inline const Stencil& getStencil(int i_r, bool DirBC_Interior) { // Only i_r = 0 is implemented. // Stencils are only used to obtain offsets to index into COO/CSR matrices. @@ -9,13 +11,39 @@ const Stencil& SmootherTake::getStencil(int i_r) const // because it has an additional across-origin coupling. assert(i_r == 0); - const bool DirBC_Interior = Smoother::DirBC_Interior_; + /* ------------------- */ + /* Stencil definitions */ + /* ------------------- */ - return DirBC_Interior ? stencil_DB_ : circle_stencil_across_origin_; + // The stencil definitions must be defined before the declaration of the inner_boundary_mumps_solver_, + // since the mumps solver will be build in the member initializer of the Smoother class. + + // Stencils encode neighborhood connectivity for A_sc matrix assembly. + // It is only used in the construction of COO/CSR matrices. + // Thus it is only used for the interior boundary matrix and not needed for the tridiagonal matrices. + // The Stencil class stores the offset for each position. + // - Non-zero matrix indicesare obtained via `ptr + offset` + // - A offset value of `-1` means the position is not included in the stencil pattern. + // - Other values (0, 1, 2, ..., stencil_size - 1) correspond to valid stencil indices. + + // clang-format off + static const Stencil stencil_DB = { + -1, -1, -1, + -1, 0, -1, + -1, -1, -1 + }; + static const Stencil circle_stencil_across_origin = { + -1, 3, -1, + 1, 0, -1, + -1, 2, -1 + }; + // clang-format on + + return DirBC_Interior ? stencil_DB : circle_stencil_across_origin; } -template -int SmootherTake::getNonZeroCountCircleAsc(int i_r) const +// Number of nonzero A_sc entries. +static inline int getNonZeroCountCircleAsc(int i_r, const PolarGrid& grid, bool DirBC_Interior) { // Only i_r = 0 is implemented. // The number of nonzero elements is only needed to construct COO matrices. @@ -23,15 +51,13 @@ int SmootherTake::getNonZeroCountCircleAsc(int i_r) const // because it has an additional across-origin coupling. assert(i_r == 0); - const PolarGrid& grid = Smoother::grid_; - const bool DirBC_Interior = Smoother::DirBC_Interior_; - const int size_stencil_inner_boundary = DirBC_Interior ? 1 : 4; return size_stencil_inner_boundary * grid.ntheta(); } -template -int SmootherTake::getCircleAscIndex(int i_r, int i_theta) const +// Obtain a ptr to index into COO matrices. +// It accumulates all stencil sizes within a line up to, but excluding the current node. +static inline int getCircleAscIndex(int i_r, int i_theta, bool DirBC_Interior) { // Only i_r = 0 is implemented. // getCircleAscIndex accumulates all stencil sizes within a line up to, but excluding the current node. @@ -40,8 +66,6 @@ int SmootherTake::getCircleAscIndex(int i_r, int i_theta) const // because it has an additional across-origin coupling. assert(i_r == 0); - const bool DirBC_Interior = Smoother::DirBC_Interior_; - const int size_stencil_inner_boundary = DirBC_Interior ? 1 : 4; return size_stencil_inner_boundary * i_theta; } diff --git a/include/Smoother/SmootherTake/smootherTake.h b/include/Smoother/SmootherTake/smootherTake.h index 63967288..1708dd81 100644 --- a/include/Smoother/SmootherTake/smootherTake.h +++ b/include/Smoother/SmootherTake/smootherTake.h @@ -62,34 +62,6 @@ class SmootherTake : public Smoother void smoothing(Vector x, ConstVector rhs, Vector temp) override; private: - /* ------------------- */ - /* Stencil definitions */ - /* ------------------- */ - - // The stencil definitions must be defined before the declaration of the inner_boundary_mumps_solver_, - // since the mumps solver will be build in the member initializer of the Smoother class. - - // Stencils encode neighborhood connectivity for A_sc matrix assembly. - // It is only used in the construction of COO/CSR matrices. - // Thus it is only used for the interior boundary matrix and not needed for the tridiagonal matrices. - // The Stencil class stores the offset for each position. - // - Non-zero matrix indicesare obtained via `ptr + offset` - // - A offset value of `-1` means the position is not included in the stencil pattern. - // - Other values (0, 1, 2, ..., stencil_size - 1) correspond to valid stencil indices. - - // clang-format off - const Stencil stencil_DB_ = { - -1, -1, -1, - -1, 0, -1, - -1, -1, -1 - }; - const Stencil circle_stencil_across_origin_ = { - -1, 3, -1, - 1, 0, -1, - -1, 2, -1 - }; - // clang-format on - /* ------------------- */ /* Tridiagonal solvers */ /* ------------------- */ @@ -128,38 +100,13 @@ class SmootherTake : public Smoother // Solver object (owns matrix if MUMPS, references if in-house solver). InnerBoundarySolver inner_boundary_solver_; - /* -------------- */ - /* Stencil access */ - /* -------------- */ - - // Select correct stencil depending on the grid position. - const Stencil& getStencil(int i_r) const; /* Only i_r = 0 implemented */ - // Number of nonzero A_sc entries. - int getNonZeroCountCircleAsc(int i_r) const; /* Only i_r = 0 implemented */ - // Obtain a ptr to index into COO matrices. - // It accumulates all stencil sizes within a line up to, but excluding the current node. - int getCircleAscIndex(int i_r, int i_theta) const; /* Only i_r = 0 implemented */ - /* --------------- */ /* Matrix assembly */ /* --------------- */ // Build all A_sc matrices for circle and radial smoothers. void buildTridiagonalSolverMatrices(); - // Build the tridiagonal solver matrices for a specific node (i_r, i_theta) - void nodeBuildTridiagonalSolverMatrices(int i_r, int i_theta, const PolarGrid& grid, bool DirBC_Interior, - BatchedTridiagonalSolver& circle_tridiagonal_solver, - BatchedTridiagonalSolver& radial_tridiagonal_solver, - ConstVector& arr, ConstVector& att, - ConstVector& art, ConstVector& detDF, - ConstVector& coeff_beta); - // Build the solver matrix for the interior boundary (i_r = 0) which is non-tridiagonal due to across-origin coupling. InnerBoundaryMatrix buildInteriorBoundarySolverMatrix(); - // Build the solver matrix for a specific node (i_r = 0, i_theta) on the interior boundary. - void nodeBuildInteriorBoundarySolverMatrix(int i_theta, const PolarGrid& grid, bool DirBC_Interior, - InnerBoundaryMatrix& matrix, ConstVector& arr, - ConstVector& att, ConstVector& art, - ConstVector& detDF, ConstVector& coeff_beta); /* ---------------------- */ /* Orthogonal application */ @@ -198,5 +145,5 @@ class SmootherTake : public Smoother #include "buildTridiagonalAsc.inl" #include "applyAscOrtho.inl" #include "solveAscSystem.inl" -#include "matrixStencil.inl" + } // namespace gmgpolar