Skip to content
Merged
Show file tree
Hide file tree
Changes from all commits
Commits
File filter

Filter by extension

Filter by extension

Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
12 changes: 0 additions & 12 deletions include/DirectSolver/DirectSolverGive/buildSolverMatrix.inl
Original file line number Diff line number Diff line change
Expand Up @@ -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;
Expand Down Expand Up @@ -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);

Expand Down Expand Up @@ -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;
Expand Down
22 changes: 9 additions & 13 deletions include/DirectSolver/DirectSolverTake/buildSolverMatrix.inl
Original file line number Diff line number Diff line change
Expand Up @@ -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 */
Expand All @@ -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;
Expand Down
13 changes: 9 additions & 4 deletions include/DirectSolver/DirectSolverTake/matrixStencil.inl
Original file line number Diff line number Diff line change
Expand Up @@ -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_ = {
Expand Down Expand Up @@ -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_;
}
}
Expand Down Expand Up @@ -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)
Expand Down Expand Up @@ -206,4 +211,4 @@ static KOKKOS_INLINE_FUNCTION bool validateSolverMatrixIndexing(const PolarGrid&
return true;
}

} // namespace direct_solver_take
} // namespace direct_solver_take
Original file line number Diff line number Diff line change
Expand Up @@ -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);
Expand Down Expand Up @@ -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);
Expand Down
Original file line number Diff line number Diff line change
Expand Up @@ -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;

Expand All @@ -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) {
Expand Down Expand Up @@ -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;

/* -------------------------- */
Expand Down
Original file line number Diff line number Diff line change
Expand Up @@ -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;
/* ------------------------------------------ */
Expand Down Expand Up @@ -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;
Expand Down Expand Up @@ -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;

Expand Down Expand Up @@ -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;

Expand Down
Original file line number Diff line number Diff line change
Expand Up @@ -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);
Expand Down
Original file line number Diff line number Diff line change
Expand Up @@ -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) {
Expand All @@ -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;
Expand All @@ -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;
Expand All @@ -172,7 +161,6 @@ ExtrapolatedSmootherTake<LevelCacheType>::buildInteriorBoundarySolverMatrix()
const LevelCacheType& level_cache = ExtrapolatedSmoother<LevelCacheType>::level_cache_;
const bool DirBC_Interior = ExtrapolatedSmoother<LevelCacheType>::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
Expand All @@ -182,6 +170,7 @@ ExtrapolatedSmootherTake<LevelCacheType>::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<double, Kokkos::HostSpace> inner_boundary_solver_matrix(ntheta, ntheta, nnz);
inner_boundary_solver_matrix.is_symmetric(true);
Expand Down
Original file line number Diff line number Diff line change
Expand Up @@ -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 */
Expand Down Expand Up @@ -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 */
Expand Down
3 changes: 1 addition & 2 deletions include/GMGPolar/solver.h
Original file line number Diff line number Diff line change
Expand Up @@ -586,8 +586,7 @@ std::pair<double, double> GMGPolar<DomainGeometry, DensityProfileCoefficients>::
Level<DomainGeometry, DensityProfileCoefficients>& level, HostConstVector<double> solution,
HostVector<double> error, const ExactSolution& exact_solution)
{
const PolarGrid& grid = level.grid();
const LevelCache<DomainGeometry, DensityProfileCoefficients>& levelCache = level.levelCache();
const PolarGrid& grid = level.grid();

assert(solution.size() == error.size());
assert(std::ssize(solution) == grid.numberOfNodes());
Expand Down
3 changes: 1 addition & 2 deletions include/GMGPolar/utils.h
Original file line number Diff line number Diff line change
Expand Up @@ -243,8 +243,7 @@ void GMGPolar<DomainGeometry, DensityProfileCoefficients>::writeToVTK(
const std::filesystem::path& file_path, const Level<DomainGeometry, DensityProfileCoefficients>& level,
HostConstVector<double> grid_function)
{
const PolarGrid& grid = level.grid();
const LevelCache<DomainGeometry, DensityProfileCoefficients>& level_cache = level.levelCache();
const PolarGrid& grid = level.grid();

assert(std::ssize(grid_function) == grid.numberOfNodes());

Expand Down
8 changes: 2 additions & 6 deletions include/Smoother/SmootherGive/buildInnerBoundaryAsc.inl
Original file line number Diff line number Diff line change
Expand Up @@ -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;

Expand Down Expand Up @@ -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) */
Expand Down Expand Up @@ -265,7 +261,6 @@ SmootherGive<LevelCacheType>::buildInteriorBoundarySolverMatrix()
const LevelCacheType& level_cache = Smoother<LevelCacheType>::level_cache_;
const bool DirBC_Interior = Smoother<LevelCacheType>::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
Expand All @@ -275,6 +270,7 @@ SmootherGive<LevelCacheType>::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<double, Kokkos::HostSpace> inner_boundary_solver_matrix(ntheta, ntheta, nnz);
inner_boundary_solver_matrix.is_symmetric(true);
Expand Down
Loading
Loading