diff --git a/include/DirectSolver/DirectSolverGive/buildSolverMatrix.inl b/include/DirectSolver/DirectSolverGive/buildSolverMatrix.inl index 796f83ec..7c2804ce 100644 --- a/include/DirectSolver/DirectSolverGive/buildSolverMatrix.inl +++ b/include/DirectSolver/DirectSolverGive/buildSolverMatrix.inl @@ -227,16 +227,11 @@ nodeBuildSolverMatrixGive(const int i_r, const int i_theta, const PolarGrid& gri 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, grid, DirBC_Interior); const int right_nz_index = getSolverMatrixIndex(i_r + 1, i_theta, grid, DirBC_Interior); 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; @@ -632,13 +627,11 @@ nodeBuildSolverMatrixGive(const int i_r, const int i_theta, const PolarGrid& gri const int center_nz_index = getSolverMatrixIndex(i_r, i_theta, grid, DirBC_Interior); const int left_nz_index = getSolverMatrixIndex(i_r - 1, i_theta, grid, DirBC_Interior); - const int right_nz_index = getSolverMatrixIndex(i_r + 1, i_theta, grid, DirBC_Interior); const int bottom_nz_index = getSolverMatrixIndex(i_r, i_theta_M1, grid, DirBC_Interior); const int top_nz_index = getSolverMatrixIndex(i_r, i_theta_P1, grid, DirBC_Interior); 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); @@ -760,16 +753,11 @@ nodeBuildSolverMatrixGive(const int i_r, const int i_theta, const PolarGrid& gri const 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, grid, DirBC_Interior); const int left_nz_index = getSolverMatrixIndex(i_r - 1, i_theta, grid, DirBC_Interior); 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; diff --git a/include/DirectSolver/DirectSolverTake/buildSolverMatrix.inl b/include/DirectSolver/DirectSolverTake/buildSolverMatrix.inl index b1289786..f8ceedc1 100644 --- a/include/DirectSolver/DirectSolverTake/buildSolverMatrix.inl +++ b/include/DirectSolver/DirectSolverTake/buildSolverMatrix.inl @@ -382,15 +382,13 @@ nodeBuildSolverMatrixTake(const int i_r, const int i_theta, const PolarGrid& gri const int center_nz_index = getSolverMatrixIndex(i_r, i_theta, grid, DirBC_Interior); - 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 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 top_left_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 */ @@ -405,10 +403,8 @@ nodeBuildSolverMatrixTake(const int i_r, const int i_theta, const PolarGrid& gri - 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 */ + const double bottom_left_value = -0.25 * (art(left_index) + art(bottom_index)); /* Bottom Left */ + const double top_left_value = +0.25 * (art(left_index) + art(top_index)); /* Top Left */ /* Fill matrix row of (i,j) */ row = center_index; diff --git a/include/DirectSolver/DirectSolverTake/matrixStencil.inl b/include/DirectSolver/DirectSolverTake/matrixStencil.inl index 6d223bb8..71df1224 100644 --- a/include/DirectSolver/DirectSolverTake/matrixStencil.inl +++ b/include/DirectSolver/DirectSolverTake/matrixStencil.inl @@ -29,12 +29,14 @@ static KOKKOS_INLINE_FUNCTION int getStencilSize(int global_index, const PolarGr else if (i_r == grid.nr() - 2) { return size_stencil_next_outer_boundary; } + Kokkos::abort("Invalid stencil index"); + return -1; } static KOKKOS_INLINE_FUNCTION const Stencil& getStencil(const int i_r, const PolarGrid& grid, const bool DirBC_Interior) { - assert(0 <= i_r && i_r < grid.nr()); - assert(grid.nr() >= 4); + KOKKOS_ASSERT(0 <= i_r && i_r < grid.nr()); + KOKKOS_ASSERT(grid.nr() >= 4); // clang-format off static constexpr Stencil stencil_interior_ = { @@ -76,7 +78,8 @@ static KOKKOS_INLINE_FUNCTION const Stencil& getStencil(const int i_r, const Pol else if (i_r == 1 && DirBC_Interior) { return stencil_next_inner_DB_; } - else if (i_r == grid.nr() - 2) { + else { + //if (i_r == grid.nr() - 2) return stencil_next_outer_DB_; } } @@ -175,6 +178,8 @@ static KOKKOS_INLINE_FUNCTION int getSolverMatrixIndex(const int i_r, const int size_stencil_next_outer_boundary * prior_next_outer_boundary_nodes + size_stencil_outer_boundary * prior_outer_boundary_nodes; } + Kokkos::abort("Invalid stencil index"); + return -1; } static KOKKOS_INLINE_FUNCTION bool validateSolverMatrixIndexing(const PolarGrid& grid, const bool DirBC_Interior) @@ -206,4 +211,4 @@ static KOKKOS_INLINE_FUNCTION bool validateSolverMatrixIndexing(const PolarGrid& return true; } -} // namespace direct_solver_take \ No newline at end of file +} // namespace direct_solver_take diff --git a/include/ExtrapolatedSmoother/ExtrapolatedSmootherGive/applyAscOrtho.inl b/include/ExtrapolatedSmoother/ExtrapolatedSmootherGive/applyAscOrtho.inl index 0ab19ccf..e1edf5cb 100644 --- a/include/ExtrapolatedSmoother/ExtrapolatedSmootherGive/applyAscOrtho.inl +++ b/include/ExtrapolatedSmoother/ExtrapolatedSmootherGive/applyAscOrtho.inl @@ -149,11 +149,11 @@ nodeApplyAscOrthoCircleGiveInside(const int i_r, const int i_theta, const PolarG 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 i_theta_Across = grid.wrapThetaIndex(i_theta + grid.ntheta() / 2); + const int i_theta_M1 = grid.wrapThetaIndex(i_theta - 1); + const int i_theta_P1 = grid.wrapThetaIndex(i_theta + 1); + //const int i_theta_Across = grid.wrapThetaIndex(i_theta + grid.ntheta() / 2); - const int left = grid.index(i_r, i_theta_Across); + //const int left = grid.index(i_r, i_theta_Across); const int bottom = grid.index(i_r, i_theta_M1); const int top = grid.index(i_r, i_theta_P1); const int right = grid.index(i_r + 1, i_theta); @@ -471,7 +471,6 @@ nodeApplyAscOrthoRadialGiveInside(const int i_r, const int i_theta, const PolarG } } else if (i_r == grid.numberSmootherCircles() - 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); diff --git a/include/ExtrapolatedSmoother/ExtrapolatedSmootherGive/buildInnerBoundaryAsc.inl b/include/ExtrapolatedSmoother/ExtrapolatedSmootherGive/buildInnerBoundaryAsc.inl index c896eabd..5ab6ce4d 100644 --- a/include/ExtrapolatedSmoother/ExtrapolatedSmootherGive/buildInnerBoundaryAsc.inl +++ b/include/ExtrapolatedSmoother/ExtrapolatedSmootherGive/buildInnerBoundaryAsc.inl @@ -94,7 +94,6 @@ nodeBuildInteriorBoundarySolverMatrix_i_r_0(const int i_theta, const PolarGrid& const int center_index = i_theta; const int left_index = i_theta_AcrossOrigin; - const int right_index = i_theta; const int bottom_index = i_theta_M1; const int top_index = i_theta_P1; @@ -103,7 +102,6 @@ nodeBuildInteriorBoundarySolverMatrix_i_r_0(const int i_theta, const PolarGrid& const int top_nz_index = getCircleAscIndex(i_r, i_theta_P1, grid, DirBC_Interior); const int left_nz_index = getCircleAscIndex(i_r, i_theta_AcrossOrigin, grid, DirBC_Interior); - int nz_index; const Stencil& CenterStencil = getStencil(i_r, i_theta, grid, DirBC_Interior); if (i_theta & 1) { @@ -219,15 +217,11 @@ nodeBuildInteriorBoundarySolverMatrix_i_r_1(const int i_theta, const PolarGrid& level_cache.obtainValues(i_r, i_theta, center, radius, theta, coeff_beta, arr, att, art, detDF); 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 int i_theta_M1 = grid.wrapThetaIndex(i_theta - 1); - const int i_theta_P1 = grid.wrapThetaIndex(i_theta + 1); - const int left_index = i_theta; /* -------------------------- */ diff --git a/include/ExtrapolatedSmoother/ExtrapolatedSmootherGive/buildTridiagonalAsc.inl b/include/ExtrapolatedSmoother/ExtrapolatedSmootherGive/buildTridiagonalAsc.inl index a2d32ebf..b044dc4a 100644 --- a/include/ExtrapolatedSmoother/ExtrapolatedSmootherGive/buildTridiagonalAsc.inl +++ b/include/ExtrapolatedSmoother/ExtrapolatedSmootherGive/buildTridiagonalAsc.inl @@ -44,7 +44,6 @@ nodeBuildTridiagonalSolverMatricesCircleSection(const int i_r, const int i_theta double coeff_beta, arr, att, art, detDF; level_cache.obtainValues(i_r, i_theta, center, radius, theta, coeff_beta, arr, att, art, detDF); - int ptr, offset; int row, column; double value; /* ------------------------------------------ */ @@ -241,13 +240,8 @@ nodeBuildTridiagonalSolverMatricesCircleSection(const int i_r, const int 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 = i_theta; const int right_index = i_theta; - const int bottom_index = i_theta_M1; - const int top_index = i_theta_P1; /* Fill matrix row of (i,j) */ row = center_index; @@ -511,7 +505,6 @@ nodeBuildTridiagonalSolverMatricesRadialSection(const int i_r, const int i_theta double coeff_beta, arr, att, art, detDF; level_cache.obtainValues(i_r, i_theta, center, radius, theta, coeff_beta, arr, att, art, detDF); - int ptr, offset; int row, column; double value; @@ -918,7 +911,7 @@ nodeBuildTridiagonalSolverMatricesRadialSection(const int i_r, const int i_theta const int center_index = i_r - numberSmootherCircles; const int left_index = i_r - numberSmootherCircles - 1; - const int right_index = i_r - numberSmootherCircles + 1; + //const int right_index = i_r - numberSmootherCircles + 1; const int bottom_index = i_r - numberSmootherCircles; const int top_index = i_r - numberSmootherCircles; diff --git a/include/ExtrapolatedSmoother/ExtrapolatedSmootherTake/applyAscOrtho.inl b/include/ExtrapolatedSmoother/ExtrapolatedSmootherTake/applyAscOrtho.inl index 96042177..97be93b6 100644 --- a/include/ExtrapolatedSmoother/ExtrapolatedSmootherTake/applyAscOrtho.inl +++ b/include/ExtrapolatedSmoother/ExtrapolatedSmootherTake/applyAscOrtho.inl @@ -135,16 +135,15 @@ nodeApplyAscOrthoCircleTake(const int i_r, const int i_theta, const PolarGrid& g 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 i_theta_Across = grid.wrapThetaIndex(i_theta + grid.ntheta() / 2); + const int i_theta_M1 = grid.wrapThetaIndex(i_theta - 1); + const int i_theta_P1 = grid.wrapThetaIndex(i_theta + 1); + //const int i_theta_Across = grid.wrapThetaIndex(i_theta + grid.ntheta() / 2); - const int left = grid.index(i_r, i_theta_Across); + //const int left = grid.index(i_r, i_theta_Across); const int bottom = grid.index(i_r, i_theta_M1); const int center = grid.index(i_r, i_theta); const int top = grid.index(i_r, i_theta_P1); diff --git a/include/ExtrapolatedSmoother/ExtrapolatedSmootherTake/buildInnerBoundaryAsc.inl b/include/ExtrapolatedSmoother/ExtrapolatedSmootherTake/buildInnerBoundaryAsc.inl index 0a5dd32c..55c71ba9 100644 --- a/include/ExtrapolatedSmoother/ExtrapolatedSmootherTake/buildInnerBoundaryAsc.inl +++ b/include/ExtrapolatedSmoother/ExtrapolatedSmootherTake/buildInnerBoundaryAsc.inl @@ -90,16 +90,9 @@ nodeBuildInteriorBoundarySolverMatrix(const int i_theta, const PolarGrid& grid, const int center_index = i_theta; const int left_index = i_theta_AcrossOrigin; - const int right_index = i_theta; - const int bottom_index = i_theta_M1; - const int top_index = i_theta_P1; const int center_nz_index = getCircleAscIndex(i_r, i_theta, grid, DirBC_Interior); - const int bottom_nz_index = getCircleAscIndex(i_r, i_theta_M1, grid, DirBC_Interior); - const int top_nz_index = getCircleAscIndex(i_r, i_theta_P1, grid, DirBC_Interior); - const int left_nz_index = getCircleAscIndex(i_r, i_theta_AcrossOrigin, grid, DirBC_Interior); - int nz_index; const Stencil& CenterStencil = getStencil(i_r, i_theta, grid, DirBC_Interior); if (i_theta & 1) { @@ -125,8 +118,6 @@ nodeBuildInteriorBoundarySolverMatrix(const int i_theta, const PolarGrid& grid, row = center_index; ptr = center_nz_index; - const Stencil& CenterStencil = getStencil(i_r, i_theta, grid, DirBC_Interior); - offset = CenterStencil[StencilPosition::Center]; column = center_index; value = center_value; @@ -149,8 +140,6 @@ nodeBuildInteriorBoundarySolverMatrix(const int i_theta, const PolarGrid& grid, row = center_index; ptr = center_nz_index; - const Stencil& CenterStencil = getStencil(i_r, i_theta, grid, DirBC_Interior); - offset = CenterStencil[StencilPosition::Center]; column = center_index; value = 1.0; @@ -172,7 +161,6 @@ ExtrapolatedSmootherTake::buildInteriorBoundarySolverMatrix() const LevelCacheType& level_cache = ExtrapolatedSmoother::level_cache_; const bool DirBC_Interior = ExtrapolatedSmoother::DirBC_Interior_; - const int i_r = 0; const int ntheta = grid.ntheta(); // The interior boundary matrix is symmetric due to the periodicity in the theta direction @@ -182,6 +170,7 @@ ExtrapolatedSmootherTake::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 i_r = 0; const int nnz = getNonZeroCountCircleAsc(i_r, grid, DirBC_Interior); SparseMatrixCOO inner_boundary_solver_matrix(ntheta, ntheta, nnz); inner_boundary_solver_matrix.is_symmetric(true); diff --git a/include/ExtrapolatedSmoother/ExtrapolatedSmootherTake/buildTridiagonalAsc.inl b/include/ExtrapolatedSmoother/ExtrapolatedSmootherTake/buildTridiagonalAsc.inl index 814d88e5..38ae34a5 100644 --- a/include/ExtrapolatedSmoother/ExtrapolatedSmootherTake/buildTridiagonalAsc.inl +++ b/include/ExtrapolatedSmoother/ExtrapolatedSmootherTake/buildTridiagonalAsc.inl @@ -32,9 +32,8 @@ static KOKKOS_INLINE_FUNCTION void nodeBuildTridiagonalSolverMatricesCircleSecti KOKKOS_ASSERT(numberSmootherCircles >= 3); KOKKOS_ASSERT(lengthRadialSmoother >= 3); - int ptr, offset; - int row, column, col; - double value, val; + int row, column; + double value; /* ------------------------------------------ */ /* Node in the interior of the Circle Section */ @@ -182,9 +181,8 @@ static KOKKOS_INLINE_FUNCTION void nodeBuildTridiagonalSolverMatricesRadialSecti KOKKOS_ASSERT(numberSmootherCircles >= 3); KOKKOS_ASSERT(lengthRadialSmoother >= 3); - int ptr, offset; - int row, column, col; - double value, val; + int row, column; + double value; /* ------------------------------------------ */ /* Node in the interior of the Radial Section */ diff --git a/include/GMGPolar/solver.h b/include/GMGPolar/solver.h index 7a2477a4..6a38ee74 100644 --- a/include/GMGPolar/solver.h +++ b/include/GMGPolar/solver.h @@ -586,8 +586,7 @@ std::pair GMGPolar:: Level& level, HostConstVector solution, HostVector error, const ExactSolution& exact_solution) { - const PolarGrid& grid = level.grid(); - const LevelCache& levelCache = level.levelCache(); + const PolarGrid& grid = level.grid(); assert(solution.size() == error.size()); assert(std::ssize(solution) == grid.numberOfNodes()); diff --git a/include/GMGPolar/utils.h b/include/GMGPolar/utils.h index ca637fa1..84646486 100644 --- a/include/GMGPolar/utils.h +++ b/include/GMGPolar/utils.h @@ -243,8 +243,7 @@ void GMGPolar::writeToVTK( const std::filesystem::path& file_path, const Level& level, HostConstVector grid_function) { - const PolarGrid& grid = level.grid(); - const LevelCache& level_cache = level.levelCache(); + const PolarGrid& grid = level.grid(); assert(std::ssize(grid_function) == grid.numberOfNodes()); diff --git a/include/Smoother/SmootherGive/buildInnerBoundaryAsc.inl b/include/Smoother/SmootherGive/buildInnerBoundaryAsc.inl index 8350bc17..4e867f4d 100644 --- a/include/Smoother/SmootherGive/buildInnerBoundaryAsc.inl +++ b/include/Smoother/SmootherGive/buildInnerBoundaryAsc.inl @@ -108,7 +108,7 @@ nodeBuildInteriorBoundarySolverMatrix_i_r_0(const int i_theta, const PolarGrid& const int top_nz_index = getCircleAscIndex(i_r, i_theta_P1, DirBC_Interior); const int left_nz_index = getCircleAscIndex(i_r, i_theta_AcrossOrigin, DirBC_Interior); - int nz_index; /* Fill matrix row of (i,j) */ + /* Fill matrix row of (i,j) */ row = center_index; ptr = center_nz_index; @@ -226,15 +226,11 @@ nodeBuildInteriorBoundarySolverMatrix_i_r_1(const int i_theta, const PolarGrid& level_cache.obtainValues(i_r, i_theta, center, radius, theta, coeff_beta, arr, att, art, detDF); 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 int i_theta_M1 = grid.wrapThetaIndex(i_theta - 1); - const int i_theta_P1 = grid.wrapThetaIndex(i_theta + 1); - const int left_index = i_theta; /* Fill matrix row of (i-1,j) */ @@ -265,7 +261,6 @@ SmootherGive::buildInteriorBoundarySolverMatrix() const LevelCacheType& level_cache = Smoother::level_cache_; const bool DirBC_Interior = Smoother::DirBC_Interior_; - const int i_r = 0; const int ntheta = grid.ntheta(); // The interior boundary matrix is symmetric due to the periodicity in the theta direction @@ -275,6 +270,7 @@ SmootherGive::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 i_r = 0; const int nnz = getNonZeroCountCircleAsc(i_r, grid, DirBC_Interior); SparseMatrixCOO inner_boundary_solver_matrix(ntheta, ntheta, nnz); inner_boundary_solver_matrix.is_symmetric(true); diff --git a/include/Smoother/SmootherGive/buildTridiagonalAsc.inl b/include/Smoother/SmootherGive/buildTridiagonalAsc.inl index 20270402..55071c6b 100644 --- a/include/Smoother/SmootherGive/buildTridiagonalAsc.inl +++ b/include/Smoother/SmootherGive/buildTridiagonalAsc.inl @@ -42,7 +42,6 @@ nodeBuildTridiagonalSolverMatricesCircleSection(const int i_r, const int i_theta double coeff_beta, arr, att, art, detDF; level_cache.obtainValues(i_r, i_theta, center, radius, theta, coeff_beta, arr, att, art, detDF); - int ptr, offset; int row, column; double value; @@ -163,9 +162,6 @@ nodeBuildTridiagonalSolverMatricesCircleSection(const int i_r, const int 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 auto& center_solver = circle_tridiagonal_solver; const int center_batch = i_r; const auto& right_solver = circle_tridiagonal_solver; @@ -173,8 +169,6 @@ nodeBuildTridiagonalSolverMatricesCircleSection(const int i_r, const int i_theta const int center_index = i_theta; const int right_index = i_theta; - const int bottom_index = i_theta_M1; - const int top_index = i_theta_P1; /* Fill matrix row of (i,j) */ row = center_index; @@ -218,7 +212,6 @@ nodeBuildTridiagonalSolverMatricesRadialSection(const int i_r, const int i_theta double coeff_beta, arr, att, art, detDF; level_cache.obtainValues(i_r, i_theta, center, radius, theta, coeff_beta, arr, att, art, detDF); - int ptr, offset; int row, column; double value; @@ -432,7 +425,6 @@ nodeBuildTridiagonalSolverMatricesRadialSection(const int i_r, const int i_theta const int center_index = i_r - numberSmootherCircles; const int left_index = i_r - numberSmootherCircles - 1; - const int right_index = i_r - numberSmootherCircles + 1; const int bottom_index = i_r - numberSmootherCircles; const int top_index = i_r - numberSmootherCircles; diff --git a/include/Smoother/SmootherTake/buildInnerBoundaryAsc.inl b/include/Smoother/SmootherTake/buildInnerBoundaryAsc.inl index 4ae3651b..d0ef443d 100644 --- a/include/Smoother/SmootherTake/buildInnerBoundaryAsc.inl +++ b/include/Smoother/SmootherTake/buildInnerBoundaryAsc.inl @@ -150,7 +150,6 @@ SmootherTake::buildInteriorBoundarySolverMatrix() const LevelCacheType& level_cache = Smoother::level_cache_; const bool DirBC_Interior = Smoother::DirBC_Interior_; - const int i_r = 0; const int ntheta = grid.ntheta(); // The interior boundary matrix is symmetric due to the periodicity in the theta direction @@ -160,6 +159,7 @@ 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 i_r = 0; const int nnz = getNonZeroCountCircleAsc(i_r, grid, DirBC_Interior); SparseMatrixCOO inner_boundary_solver_matrix(ntheta, ntheta, nnz); inner_boundary_solver_matrix.is_symmetric(true); diff --git a/src/InputFunctions/SourceTerms/cartesianR2_Poisson_CzarnyGeometry.cpp b/src/InputFunctions/SourceTerms/cartesianR2_Poisson_CzarnyGeometry.cpp index 68056746..f13acfbf 100644 --- a/src/InputFunctions/SourceTerms/cartesianR2_Poisson_CzarnyGeometry.cpp +++ b/src/InputFunctions/SourceTerms/cartesianR2_Poisson_CzarnyGeometry.cpp @@ -26,7 +26,6 @@ KOKKOS_FUNCTION double CartesianR2_Poisson_CzarnyGeometry::operator()(std::size_ double temp = sqrt(inverse_aspect_ratio_epsilon * (inverse_aspect_ratio_epsilon + 2.0 * (r / Rmax) * cos_theta) + 1.0); double sin_theta_pow2 = pow(sin_theta, 2.0); - double cos_theta_pow2 = pow(cos_theta, 2.0); double temp1 = pow((2.0 - temp), 2.0); double temp2 = pow((2.0 - temp), 3.0); diff --git a/src/InputFunctions/SourceTerms/polarR6_ZoniShiftedGyro_CulhamGeometry.cpp b/src/InputFunctions/SourceTerms/polarR6_ZoniShiftedGyro_CulhamGeometry.cpp index 0b734bc8..54ee03b7 100644 --- a/src/InputFunctions/SourceTerms/polarR6_ZoniShiftedGyro_CulhamGeometry.cpp +++ b/src/InputFunctions/SourceTerms/polarR6_ZoniShiftedGyro_CulhamGeometry.cpp @@ -9,9 +9,7 @@ PolarR6_ZoniShiftedGyro_CulhamGeometry::PolarR6_ZoniShiftedGyro_CulhamGeometry(P KOKKOS_FUNCTION double PolarR6_ZoniShiftedGyro_CulhamGeometry::operator()(std::size_t i_r, std::size_t i_theta) const { - double r = grid_.radius(i_r); - double theta = grid_.theta(i_theta); - double sin_theta = std::sin(theta); - double cos_theta = std::cos(theta); + double r = grid_.radius(i_r); + double theta = grid_.theta(i_theta); return 0.4096 * pow((r / Rmax), 6.0) * pow(((r / Rmax) - 1.0), 6.0) * cos(11.0 * theta); -} \ No newline at end of file +} diff --git a/src/InputFunctions/SourceTerms/refined_ZoniShiftedGyro_CulhamGeometry.cpp b/src/InputFunctions/SourceTerms/refined_ZoniShiftedGyro_CulhamGeometry.cpp index 41933510..2d0b4ee1 100644 --- a/src/InputFunctions/SourceTerms/refined_ZoniShiftedGyro_CulhamGeometry.cpp +++ b/src/InputFunctions/SourceTerms/refined_ZoniShiftedGyro_CulhamGeometry.cpp @@ -9,14 +9,12 @@ Refined_ZoniShiftedGyro_CulhamGeometry::Refined_ZoniShiftedGyro_CulhamGeometry(P KOKKOS_FUNCTION double Refined_ZoniShiftedGyro_CulhamGeometry::operator()(std::size_t i_r, std::size_t i_theta) const { - double r = grid_.radius(i_r); - double theta = grid_.theta(i_theta); - double sin_theta = std::sin(theta); - double cos_theta = std::cos(theta); + double r = grid_.radius(i_r); + double theta = grid_.theta(i_theta); return ((-3.33823779536505e-15) * ((r / Rmax) * (r / Rmax)) - 0.0 * (r / Rmax) - 0.0 + exp((-3333.33333333333) * pow(((r / Rmax) - 0.9), 2.0))) * cos(21.0 * theta) + (0.00184273372222541 * ((r / Rmax) * (r / Rmax)) - 0.0018029383826828 * (r / Rmax) - 4.00652973929511e-05 + exp((-50.0) * pow(((r / Rmax) - 0.45), 2.0))) * cos(9.0 * theta); -} \ No newline at end of file +} diff --git a/src/main.cpp b/src/main.cpp index e02a35e6..e8e535f8 100644 --- a/src/main.cpp +++ b/src/main.cpp @@ -72,7 +72,6 @@ int main(int argc, char* argv[]) // --- Retrieve solution and associated grid --- // HostVector solution = solver.solution(); - const PolarGrid& grid = solver.grid(); // Print timing statistics for each solver phase solver.printTimings();