From ef2a35c394087bf63ca370c001d46774114ac983 Mon Sep 17 00:00:00 2001 From: julianlitz Date: Mon, 17 Aug 2026 00:14:46 +0200 Subject: [PATCH 1/7] Faster Take kernels --- include/Definitions/geometry_helper.h | 3 + .../DirectSolverGive/applySymmetryShift.inl | 16 +- .../DirectSolverGive/buildSolverMatrix.inl | 48 +-- .../DirectSolverTake/applySymmetryShift.inl | 8 +- .../DirectSolverTake/buildSolverMatrix.inl | 24 +- .../applyAscOrtho.inl | 84 ++--- .../applyAscOrtho.inl | 72 ++-- include/Residual/ResidualGive/applyAGive.inl | 40 +-- include/Residual/ResidualTake/applyATake.inl | 140 +++++--- .../Smoother/SmootherGive/applyAscOrtho.inl | 32 +- .../Smoother/SmootherTake/applyAscOrtho.inl | 328 +++++++++++++----- 11 files changed, 512 insertions(+), 283 deletions(-) diff --git a/include/Definitions/geometry_helper.h b/include/Definitions/geometry_helper.h index 1245e598..5f04987b 100644 --- a/include/Definitions/geometry_helper.h +++ b/include/Definitions/geometry_helper.h @@ -33,4 +33,7 @@ KOKKOS_INLINE_FUNCTION void compute_jacobian_elements(const DomainGeometry& doma /* 1.0 / det(DF) * */ /* [Jtt, -Jrt] */ /* [-Jtr, Jrr] */ + /* The stencil coefficient art is scaled by 0.25 to account for the + /* 9-point stencil used in the finite difference approximation */ + art *= 0.25; } diff --git a/include/DirectSolver/DirectSolverGive/applySymmetryShift.inl b/include/DirectSolver/DirectSolverGive/applySymmetryShift.inl index 56f48a9e..13aa872d 100644 --- a/include/DirectSolver/DirectSolverGive/applySymmetryShift.inl +++ b/include/DirectSolver/DirectSolverGive/applySymmetryShift.inl @@ -45,8 +45,8 @@ void DirectSolverGive::applySymmetryShiftInnerBoundary(Vector::applySymmetryShiftInnerBoundary(Vector::applySymmetryShiftOuterBoundary(Vector::applySymmetryShiftOuterBoundary(Vector::applySymmetryShiftInnerBoundary(Vector::applySymmetryShiftOuterBoundary(Vector 1 || (i_r == 1 && !DirBC_Interior)) { /* Don't give to the inner dirichlet boundary! */ result[left] += (-coeff1 * arr * x[center] /* Right */ + coeff1 * arr * x[left] /* Center: (Right) */ - - 0.25 * art * x[top] /* Top Right */ - + 0.25 * art * x[bottom]); /* Bottom Right */ + - art * x[top] /* Top Right */ + + art * x[bottom]); /* Bottom Right */ } /* Fill result(i+1,j) */ if (i_r < grid.nr() - 2) { /* Don't give to the outer dirichlet boundary! */ result[right] += (-coeff2 * arr * x[center] /* Left */ + coeff2 * arr * x[right] /* Center: (Left) */ - + 0.25 * art * x[top] /* Top Left */ - - 0.25 * art * x[bottom]); /* Bottom Left */ + + art * x[top] /* Top Left */ + - art * x[bottom]); /* Bottom Left */ } /* Fill result(i,j-1) */ result[bottom] += (-coeff3 * att * x[center] /* Top */ + coeff3 * att * x[bottom] /* Center: (Top) */ - - 0.25 * art * x[right] /* Top Right */ - + 0.25 * art * x[left]); /* Top Left */ + - art * x[right] /* Top Right */ + + art * x[left]); /* Top Left */ /* Fill result(i,j+1) */ result[top] += (-coeff4 * att * x[center] /* Bottom */ + coeff4 * att * x[top] /* Center: (Bottom) */ - + 0.25 * art * x[right] /* Bottom Right */ - - 0.25 * art * x[left]); /* Bottom Left */ + + art * x[right] /* Bottom Right */ + - art * x[left]); /* Bottom Left */ } /* -------------------------- */ /* Node on the inner boundary */ @@ -102,8 +102,8 @@ static KOKKOS_INLINE_FUNCTION void node_apply_a_give(int i_r, int i_theta, const /* Fill result(i+1,j) */ result[right] += (-coeff2 * arr * x[center] /* Left */ + coeff2 * arr * x[right] /* Center: (Left) */ - + 0.25 * art * x[top] /* Top Left */ - - 0.25 * art * x[bottom]); /* Bottom Left */ + + art * x[top] /* Top Left */ + - art * x[bottom]); /* Bottom Left */ } else { /* ------------------------------------------------------------- */ @@ -144,23 +144,23 @@ static KOKKOS_INLINE_FUNCTION void node_apply_a_give(int i_r, int i_theta, const /* From view the view of the across origin node, the directions are roatated by 180 degrees in the stencil! */ result[left] += (-coeff1 * arr * x[center] /* Right -> Left */ + coeff1 * arr * x[left]); /* Center: (Right) -> Center: (Left)*/ - /* + 0.25 * art * x[top]; // Top Right -> Bottom Left: REMOVED DUE TO ARTIFICAL 7 POINT STENCIL */ - /* - 0.25 * art * x[bottom]; // Bottom Right -> Top Left: REMOVED DUE TO ARTIFICAL 7 POINT STENCIL */ + /* + art * x[top]; // Top Right -> Bottom Left: REMOVED DUE TO ARTIFICAL 7 POINT STENCIL */ + /* - art * x[bottom]; // Bottom Right -> Top Left: REMOVED DUE TO ARTIFICAL 7 POINT STENCIL */ /* Fill result(i+1,j) */ result[right] += (-coeff2 * arr * x[center] /* Left */ + coeff2 * arr * x[right] /* Center: (Left) */ - + 0.25 * art * x[top] /* Top Left */ - - 0.25 * art * x[bottom]); /* Bottom Left */ + + art * x[top] /* Top Left */ + - art * x[bottom]); /* Bottom Left */ /* Fill result(i,j-1) */ result[bottom] += (-coeff3 * att * x[center] /* Top */ + coeff3 * att * x[bottom] /* Center: (Top) */ - - 0.25 * art * x[right]); /* Top Right */ - /* + 0.25 * art * x[left]; // Top Left: REMOVED DUE TO ARTIFICAL 7 POINT STENCIL */ + - art * x[right]); /* Top Right */ + /* + art * x[left]; // Top Left: REMOVED DUE TO ARTIFICAL 7 POINT STENCIL */ /* Fill result(i,j+1) */ result[top] += (-coeff4 * att * x[center] /* Bottom */ + coeff4 * att * x[top] /* Center: (Bottom) */ - + 0.25 * art * x[right]); /* Bottom Right */ - /* - 0.25 * art * x[left]; // Bottom Left: REMOVED DUE TO ARTIFICAL 7 POINT STENCIL */ + + art * x[right]); /* Bottom Right */ + /* - art * x[left]; // Bottom Left: REMOVED DUE TO ARTIFICAL 7 POINT STENCIL */ } } /* ----------------------------- */ @@ -187,8 +187,8 @@ static KOKKOS_INLINE_FUNCTION void node_apply_a_give(int i_r, int i_theta, const /* Fill result(i-1,j) */ result[left] += (-coeff1 * arr * x[center] /* Right */ + coeff1 * arr * x[left] /* Center: (Right) */ - - 0.25 * art * x[top] /* Top Right */ - + 0.25 * art * x[bottom]); /* Bottom Right */ + - art * x[top] /* Top Right */ + + art * x[bottom]); /* Bottom Right */ } } diff --git a/include/Residual/ResidualTake/applyATake.inl b/include/Residual/ResidualTake/applyATake.inl index 9d39f520..daeb5b12 100644 --- a/include/Residual/ResidualTake/applyATake.inl +++ b/include/Residual/ResidualTake/applyATake.inl @@ -3,23 +3,31 @@ namespace residual_take { -static KOKKOS_INLINE_FUNCTION void node_apply_a_take(const int i_r, const int i_theta, const PolarGrid& grid, - bool DirBC_Interior, Vector& result, - ConstVector& x, ConstVector& arr, - ConstVector& att, ConstVector& art, - ConstVector& detDF, ConstVector& coeff_beta) +static KOKKOS_INLINE_FUNCTION void applySystemMatrixTakeInterior(const int i_r, const int i_theta, + const PolarGrid& grid, bool DirBC_Interior, + Vector& result, ConstVector& x, + ConstVector& arr, ConstVector& att, + ConstVector& art, ConstVector& detDF, + ConstVector& coeff_beta) { - const int center = grid.index(i_r, i_theta); + KOKKOS_ASSERT(0 < i_r && i_r < grid.nr() - 1); - if ((i_r == 0 && DirBC_Interior) || (i_r == grid.nr() - 1)) { - result[center] = x[center]; - return; - } + const int i_theta_M1 = grid.wrapThetaIndex(i_theta - 1); + const int i_theta_P1 = grid.wrapThetaIndex(i_theta + 1); - // Across origin: h1 gets replaced with 2 * R0. - const double h1 = (i_r == 0) ? 2.0 * grid.radius(0) : grid.radialSpacing(i_r - 1); + const int bottom_left = grid.index(i_r - 1, i_theta_M1); + const int bottom = grid.index(i_r, i_theta_M1); + const int bottom_right = grid.index(i_r + 1, i_theta_M1); + const int left = grid.index(i_r - 1, i_theta); + const int center = grid.index(i_r, i_theta); + const int right = grid.index(i_r + 1, i_theta); + const int top_left = grid.index(i_r - 1, i_theta_P1); + const int top = grid.index(i_r, i_theta_P1); + const int top_right = grid.index(i_r + 1, i_theta_P1); + + 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 k1 = grid.angularSpacing(i_theta_M1); const double k2 = grid.angularSpacing(i_theta); const double coeff1 = 0.5 * (k1 + k2) / h1; @@ -28,38 +36,67 @@ static KOKKOS_INLINE_FUNCTION void node_apply_a_take(const int i_r, const int i_ const double coeff4 = 0.5 * (h1 + h2) / k2; const double coeff5 = 0.25 * (h1 + h2) * (k1 + k2); + result[center] = (+coeff5 * coeff_beta[center] * Kokkos::fabs(detDF[center]) * x[center]; /* beta_{i,j} */ + + +coeff1 * (arr[center] + arr[left]) * (x[center] - x[left]); /* Center: (Left) - Left */ + +coeff2 * (arr[center] + arr[right]) * (x[center] - x[right]); /* Center: (Right) - Right */ + +coeff3 * (att[center] + att[bottom]) * (x[center] - x[bottom]); /* Center: (Bottom) - Bottom */ + +coeff4 * (att[center] + att[top]) * (x[center] - x[top]); /* Center: (Top) - Top */ + + -(art[left] + art[bottom]) * x[bottom_left]; /* Bottom Left */ + +(art[left] + art[top]) * x[top_left]; /* Top Left */ + +(art[right] + art[bottom]) * x[bottom_right]; /* Bottom Right */ + -(art[right] + art[top]) * x[top_right]; /* Top Right */ + ); +} + +static KOKKOS_INLINE_FUNCTION void applySystemMatrixTakeBoundary(const int i_r, const int i_theta, + const PolarGrid& grid, bool DirBC_Interior, + Vector& result, ConstVector& x, + ConstVector& arr, ConstVector& att, + ConstVector& art, ConstVector& detDF, + ConstVector& coeff_beta) +{ + KOKKOS_ASSERT(i_r == 0 || i_r == grid.nr() - 1); + + if ((i_r == 0 && DirBC_Interior) || (i_r == grid.nr() - 1)) { + const int center = grid.index(i_r, i_theta); + result[center] = x[center]; + return; + } + 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 = (i_r == 0) ? grid.index(i_r, i_theta_Across) : grid.index(i_r - 1, i_theta); - const int right = grid.index(i_r + 1, i_theta); + 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); + const int right = grid.index(i_r + 1, i_theta); - double value = 0.0; - - value += coeff5 * coeff_beta[center] * std::fabs(detDF[center]) * x[center]; /* beta_{i,j} */ + // Across origin: h1 gets replaced with 2 * R0. + const double h1 = 2.0 * grid.radius(0); + const double h2 = grid.radialSpacing(i_r); + const double k1 = grid.angularSpacing(i_theta - 1); + const double k2 = grid.angularSpacing(i_theta); - value += coeff1 * (arr[center] + arr[left]) * (x[center] - x[left]); /* Center: (Left) - Left */ - value += coeff2 * (arr[center] + arr[right]) * (x[center] - x[right]); /* Center: (Right) - Right */ - value += coeff3 * (att[center] + att[bottom]) * (x[center] - x[bottom]); /* Center: (Bottom) - Bottom */ - value += coeff4 * (att[center] + att[top]) * (x[center] - x[top]); /* Center: (Top) - Top */ + 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 double coeff5 = 0.25 * (h1 + h2) * (k1 + k2); - // Across origin: Reduce 9-point stencil to the artifical 7-point stencil. - if (i_r > 0) { - const int bottom_left = grid.index(i_r - 1, i_theta_M1); - const int top_left = grid.index(i_r - 1, i_theta_P1); - value += -0.25 * (art[left] + art[bottom]) * x[bottom_left]; /* Bottom Left */ - value += +0.25 * (art[left] + art[top]) * x[top_left]; /* Top Left */ - } + result[center] = (+coeff5 * coeff_beta[center] * Kokkos::fabs(detDF[center]) * x[center]; /* beta_{i,j} */ - const int bottom_right = grid.index(i_r + 1, i_theta_M1); - const int top_right = grid.index(i_r + 1, i_theta_P1); - value += +0.25 * (art[right] + art[bottom]) * x[bottom_right]; /* Bottom Right */ - value += -0.25 * (art[right] + art[top]) * x[top_right]; /* Top Right */ + +coeff1 * (arr[center] + arr[left]) * (x[center] - x[left]); /* Center: (Left) - Left */ + +coeff2 * (arr[center] + arr[right]) * (x[center] - x[right]); /* Center: (Right) - Right */ + +coeff3 * (att[center] + att[bottom]) * (x[center] - x[bottom]); /* Center: (Bottom) - Bottom */ + +coeff4 * (att[center] + att[top]) * (x[center] - x[top]); /* Center: (Top) - Top */ - result[center] = value; + +(art[right] + art[bottom]) * x[bottom_right]; /* Bottom Right */ + -(art[right] + art[top]) * x[top_right]; /* Top Right */ + ); } } // namespace residual_take @@ -67,7 +104,8 @@ static KOKKOS_INLINE_FUNCTION void node_apply_a_take(const int i_r, const int i_ template void ResidualTake::applySystemOperator(Vector result, ConstVector x) const { - using residual_take::node_apply_a_take; + using residual_take::applySystemMatrixTakeBoundary; + using residual_take::applySystemMatrixTakeInterior; assert(result.size() == x.size()); @@ -88,26 +126,50 @@ void ResidualTake::applySystemOperator(Vector result, Co // The For loop matches circular access pattern */ Kokkos::parallel_for( - "Residual Take: Apply System Operator (Circular)", + "Residual Take: Apply System Operator Boundary (Circular)", Kokkos::MDRangePolicy>( // Rank of the index space {0, 0}, // Starting point of the index space + {1, 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) { + applySystemMatrixTakeBoundary(i_r, i_theta, grid, DirBC_Interior, result, x, arr, att, art, detDF, + coeff_beta); + }); + Kokkos::parallel_for( + "Residual Take: Apply System Operator Interior (Circular)", + Kokkos::MDRangePolicy>( // Rank of the index space + {1, 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) { - node_apply_a_take(i_r, i_theta, grid, DirBC_Interior, result, x, arr, att, art, detDF, coeff_beta); + applySystemMatrixTakeInterior(i_r, i_theta, grid, DirBC_Interior, result, x, arr, att, art, detDF, + coeff_beta); }); /* For loop matches radial access pattern */ Kokkos::parallel_for( - "Residual Take: Apply System Operator (Radial)", + "Residual Take: Apply System Operator Interior (Radial)", Kokkos::MDRangePolicy>( // Rank of the index space {0, grid.numberSmootherCircles()}, // Starting point of the index space + {grid.ntheta(), grid.nr() - 1} // 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) { + applySystemMatrixTakeInterior(i_r, i_theta, grid, DirBC_Interior, result, x, arr, att, art, detDF, + coeff_beta); + }); + Kokkos::parallel_for( + "Residual Take: Apply System Operator Boundary (Radial)", + Kokkos::MDRangePolicy>( // Rank of the index space + {0, grid.nr() - 1}, // 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) { - node_apply_a_take(i_r, i_theta, grid, DirBC_Interior, result, x, arr, att, art, detDF, coeff_beta); + applySystemMatrixTakeBoundary(i_r, i_theta, grid, DirBC_Interior, result, x, arr, att, art, detDF, + coeff_beta); }); Kokkos::fence(); diff --git a/include/Smoother/SmootherGive/applyAscOrtho.inl b/include/Smoother/SmootherGive/applyAscOrtho.inl index 9e0b29a5..fda04684 100644 --- a/include/Smoother/SmootherGive/applyAscOrtho.inl +++ b/include/Smoother/SmootherGive/applyAscOrtho.inl @@ -40,16 +40,16 @@ static KOKKOS_INLINE_FUNCTION void nodeApplyAscOrthoCircleGiveInside(int i_r, in const double coeff1 = 0.5 * (k1 + k2) / h1; const int left = grid.index(i_r - 1, i_theta); result[center] -= (-coeff1 * arr * x[left]); /* Left */ - result[bottom] -= (+0.25 * art * x[left]); /* Top Left */ - result[top] -= (-0.25 * art * x[left]); /* Bottom Left */ + result[bottom] -= (+art * x[left]); /* Top Left */ + result[top] -= (-art * x[left]); /* Bottom Left */ } const double h2 = grid.radialSpacing(i_r); const double coeff2 = 0.5 * (k1 + k2) / h2; const int right = grid.index(i_r + 1, i_theta); result[center] -= (-coeff2 * arr * x[right]); /* Right */ - result[bottom] -= (-0.25 * art * x[right]); /* Top Right */ - result[top] -= (+0.25 * art * x[right]); /* Bottom Right */ + result[bottom] -= (-art * x[right]); /* Top Right */ + result[top] -= (+art * x[right]); /* Bottom Right */ } template @@ -85,8 +85,8 @@ static KOKKOS_INLINE_FUNCTION void nodeApplyAscOrthoCircleGiveOutside(int i_r, i const double coeff1 = 0.5 * (k1 + k2) / h1; const int left = grid.index(i_r - 1, i_theta); result[left] -= (-coeff1 * arr * x[center] /* Right */ - - 0.25 * art * x[top] /* Top Right */ - + 0.25 * art * x[bottom]); /* Bottom Right */ + - art * x[top] /* Top Right */ + + art * x[bottom]); /* Bottom Right */ } /* Fill result(i+1,j) */ if (i_r < grid.numberSmootherCircles() - 1) { @@ -94,8 +94,8 @@ static KOKKOS_INLINE_FUNCTION void nodeApplyAscOrthoCircleGiveOutside(int i_r, i const double coeff2 = 0.5 * (k1 + k2) / h2; const int right = grid.index(i_r + 1, i_theta); result[right] -= (-coeff2 * arr * x[center] /* Left */ - + 0.25 * art * x[top] /* Top Left */ - - 0.25 * art * x[bottom]); /* Bottom Left */ + + art * x[top] /* Top Left */ + - art * x[bottom]); /* Bottom Left */ } } @@ -138,14 +138,14 @@ static KOKKOS_INLINE_FUNCTION void nodeApplyAscOrthoRadialGiveInside(int i_r, in if (grid.numberSmootherCircles() <= i_r && i_r <= grid.nr() - 1) { const int left = grid.index(i_r - 1, i_theta); - result[left] -= (-0.25 * art * x[top] /* Top Right */ - + 0.25 * art * x[bottom]); /* Bottom Right */ + result[left] -= (-art * x[top] /* Top Right */ + + art * x[bottom]); /* Bottom Right */ } if (grid.numberSmootherCircles() - 1 <= i_r && i_r < grid.nr() - 2) { const int right = grid.index(i_r + 1, i_theta); - result[right] -= (+0.25 * art * x[top] /* Top Left */ - - 0.25 * art * x[bottom]); /* Bottom Left */ + result[right] -= (+art * x[top] /* Top Left */ + - art * x[bottom]); /* Bottom Left */ } if (i_r == grid.numberSmootherCircles() - 1) { @@ -214,12 +214,12 @@ static KOKKOS_INLINE_FUNCTION void nodeApplyAscOrthoRadialGiveOutside(int i_r, i /* Fill result(i,j-1) */ result[bottom] -= (-coeff3 * att * x[center] /* Top */ - - 0.25 * art * x[right] /* Top Right */ - + 0.25 * art * x[left]); /* Top Left */ + - art * x[right] /* Top Right */ + + art * x[left]); /* Top Left */ /* Fill result(i,j+1) */ result[top] -= (-coeff4 * att * x[center] /* Bottom */ - + 0.25 * art * x[right] /* Bottom Right */ - - 0.25 * art * x[left]); /* Bottom Left */ + + art * x[right] /* Bottom Right */ + - art * x[left]); /* Bottom Left */ } } diff --git a/include/Smoother/SmootherTake/applyAscOrtho.inl b/include/Smoother/SmootherTake/applyAscOrtho.inl index d8b0af86..41bcb986 100644 --- a/include/Smoother/SmootherTake/applyAscOrtho.inl +++ b/include/Smoother/SmootherTake/applyAscOrtho.inl @@ -3,101 +3,159 @@ namespace smoother_take { -static KOKKOS_INLINE_FUNCTION void nodeApplyAscOrthoCircleTake(const int i_r, const int i_theta, const PolarGrid& grid, - bool DirBC_Interior, ConstVector& x, - ConstVector& rhs, Vector& result, - ConstVector& arr, ConstVector& att, - ConstVector& art) +static KOKKOS_INLINE_FUNCTION void +nodeApplyAscOrthoCircleTakeInterior(const int i_r, const int i_theta, const PolarGrid& grid, bool DirBC_Interior, + ConstVector& x, ConstVector& rhs, Vector& result, + ConstVector& arr, ConstVector& att, ConstVector& art) { - KOKKOS_ASSERT(i_r >= 0 && i_r < grid.numberSmootherCircles()); + KOKKOS_ASSERT(i_r >= 1 && i_r < grid.numberSmootherCircles()); - const int center = grid.index(i_r, i_theta); - - result[center] = rhs[center]; + const int i_theta_M1 = grid.wrapThetaIndex(i_theta - 1); + const int i_theta_P1 = grid.wrapThetaIndex(i_theta + 1); - if (i_r == 0 && DirBC_Interior) - return; + const int bottom_left = grid.index(i_r - 1, i_theta_M1); + const int left = grid.index(i_r - 1, i_theta); + const int top_left = grid.index(i_r - 1, i_theta_P1); + 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); + const int bottom_right = grid.index(i_r + 1, i_theta_M1); + const int right = grid.index(i_r + 1, i_theta); + const int top_right = grid.index(i_r + 1, i_theta_P1); - const double k1 = grid.angularSpacing(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 int i_theta_M1 = grid.wrapThetaIndex(i_theta - 1); - const int i_theta_P1 = grid.wrapThetaIndex(i_theta + 1); + const double coeff1 = 0.5 * (k1 + k2) / h1; + const double coeff2 = 0.5 * (k1 + k2) / h2; - const int bottom = grid.index(i_r, i_theta_M1); - const int top = grid.index(i_r, i_theta_P1); + result[center] = rhs[center] - (-coeff1 * (arr[center] + arr[left]) * x[left] /* Left */ + - coeff2 * (arr[center] + arr[right]) * x[right] /* Right */ - if (i_r > 0) { - const double h1 = grid.radialSpacing(i_r - 1); - const double coeff1 = 0.5 * (k1 + k2) / h1; + - (art[left] + art[bottom]) * x[bottom_left] /* Bottom Left */ + + (art[left] + art[top]) * x[top_left] /* Top Left */ + + (art[right] + art[bottom]) * x[bottom_right] /* Bottom Right */ + - (art[right] + art[top]) * x[top_right] /* Top Right */ + ); +} - const int bottom_left = grid.index(i_r - 1, i_theta_M1); - const int left = grid.index(i_r - 1, i_theta); - const int top_left = grid.index(i_r - 1, i_theta_P1); +static KOKKOS_INLINE_FUNCTION void +nodeApplyAscOrthoCircleTakeBoundary(const int i_r, const int i_theta, const PolarGrid& grid, bool DirBC_Interior, + ConstVector& x, ConstVector& rhs, Vector& result, + ConstVector& arr, ConstVector& att, ConstVector& art) +{ + KOKKOS_ASSERT(i_r == 0); - result[center] -= (-coeff1 * (arr[center] + arr[left]) * x[left] /* Left */ - - 0.25 * (art[left] + art[bottom]) * x[bottom_left] /* Bottom Left */ - + 0.25 * (art[left] + art[top]) * x[top_left] /* Top Left */ - ); + if (DirBC_Interior) { + const int center = grid.index(i_r, i_theta); + result[center] = rhs[center]; + return; } - const double h2 = grid.radialSpacing(i_r); - 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 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); const int bottom_right = grid.index(i_r + 1, i_theta_M1); const int right = grid.index(i_r + 1, i_theta); const int top_right = grid.index(i_r + 1, i_theta_P1); - result[center] -= (-coeff2 * (arr[center] + arr[right]) * x[right] /* Right */ - + 0.25 * (art[right] + art[bottom]) * x[bottom_right] /* Bottom Right */ - - 0.25 * (art[right] + art[top]) * x[top_right] /* Top Right */ + const double h2 = grid.radialSpacing(i_r); + const double k1 = grid.angularSpacing(i_theta_M1); + const double k2 = grid.angularSpacing(i_theta); - ); + const double coeff2 = 0.5 * (k1 + k2) / h2; + + result[center] = rhs[center] - (-coeff2 * (arr[center] + arr[right]) * x[right] /* Right */ + + (art[right] + art[bottom]) * x[bottom_right] /* Bottom Right */ + - (art[right] + art[top]) * x[top_right] /* Top Right */ + ); } -static KOKKOS_INLINE_FUNCTION void nodeApplyAscOrthoRadialTake(const int i_r, const int i_theta, const PolarGrid& grid, - bool DirBC_Interior, ConstVector& x, - ConstVector& rhs, Vector& result, - ConstVector& arr, ConstVector& att, - ConstVector& art) +static KOKKOS_INLINE_FUNCTION void +nodeApplyAscOrthoRadialTakeInterior(const int i_r, const int i_theta, const PolarGrid& grid, bool DirBC_Interior, + ConstVector& x, ConstVector& rhs, Vector& result, + ConstVector& arr, ConstVector& att, ConstVector& art) { - KOKKOS_ASSERT(i_r >= grid.numberSmootherCircles() && i_r < grid.nr()); - - const int center = grid.index(i_r, i_theta); + KOKKOS_ASSERT(i_r > grid.numberSmootherCircles() && i_r < grid.nr() - 2); - result[center] = rhs[center]; + const int i_theta_M1 = grid.wrapThetaIndex(i_theta - 1); + const int i_theta_P1 = grid.wrapThetaIndex(i_theta + 1); - if (i_r == grid.nr() - 1) - return; + const int bottom_left = grid.index(i_r - 1, i_theta_M1); + const int bottom = grid.index(i_r, i_theta_M1); + const int bottom_right = grid.index(i_r + 1, i_theta_M1); + const int left = grid.index(i_r - 1, i_theta); + const int center = grid.index(i_r, i_theta); + const int right = grid.index(i_r + 1, i_theta); + const int top_left = grid.index(i_r - 1, i_theta_P1); + const int top = grid.index(i_r, i_theta_P1); + const int top_right = grid.index(i_r + 1, i_theta_P1); 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 k1 = grid.angularSpacing(i_theta_M1); const double k2 = grid.angularSpacing(i_theta); const double coeff3 = 0.5 * (h1 + h2) / k1; const double coeff4 = 0.5 * (h1 + h2) / k2; + result[center] = rhs[center] - (-coeff3 * (att[center] + att[bottom]) * x[bottom] /* Bottom */ + - coeff4 * (att[center] + att[top]) * x[top] /* Top */ + + - (art[left] + art[bottom]) * x[bottom_left] /* Bottom Left */ + + (art[right] + art[bottom]) * x[bottom_right] /* Bottom Right */ + + (art[left] + art[top]) * x[top_left] /* Top Left */ + - (art[right] + art[top]) * x[top_right] /* Top Right */ + ); +} + +static KOKKOS_INLINE_FUNCTION void +nodeApplyAscOrthoRadialTakeBoundary(const int i_r, const int i_theta, const PolarGrid& grid, bool DirBC_Interior, + ConstVector& x, ConstVector& rhs, Vector& result, + ConstVector& arr, ConstVector& att, ConstVector& art) +{ + KOKKOS_ASSERT(i_r == grid.numberSmootherCircles() || i_r == grid.nr() - 2 || i_r == grid.nr() - 1); + + if (i_r == grid.nr() - 1) { + const int center = grid.index(i_r, i_theta); + result[center] = rhs[center]; + return; + } + const int i_theta_M1 = grid.wrapThetaIndex(i_theta - 1); const int i_theta_P1 = grid.wrapThetaIndex(i_theta + 1); const int bottom_left = grid.index(i_r - 1, i_theta_M1); - const int left = grid.index(i_r - 1, i_theta); - const int top_left = grid.index(i_r - 1, i_theta_P1); const int bottom = grid.index(i_r, i_theta_M1); - const int top = grid.index(i_r, i_theta_P1); const int bottom_right = grid.index(i_r + 1, i_theta_M1); + const int left = grid.index(i_r - 1, i_theta); + const int center = grid.index(i_r, i_theta); const int right = grid.index(i_r + 1, i_theta); + const int top_left = grid.index(i_r - 1, i_theta_P1); + const int top = grid.index(i_r, i_theta_P1); const int top_right = grid.index(i_r + 1, i_theta_P1); - result[center] -= (-coeff3 * (att[center] + att[bottom]) * x[bottom] /* Bottom */ - - coeff4 * (att[center] + att[top]) * x[top] /* Top */ + 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); - - 0.25 * (art[left] + art[bottom]) * x[bottom_left] /* Bottom Left */ - + 0.25 * (art[right] + art[bottom]) * x[bottom_right] /* Bottom Right */ - + 0.25 * (art[left] + art[top]) * x[top_left] /* Top Left */ - - 0.25 * (art[right] + art[top]) * x[top_right] /* Top Right */ - ); + const double coeff3 = 0.5 * (h1 + h2) / k1; + const double coeff4 = 0.5 * (h1 + h2) / k2; + + result[center] = rhs[center] - (-coeff3 * (att[center] + att[bottom]) * x[bottom] /* Bottom */ + - coeff4 * (att[center] + att[top]) * x[top] /* Top */ + + - (art[left] + art[bottom]) * x[bottom_left] /* Bottom Left */ + + (art[right] + art[bottom]) * x[bottom_right] /* Bottom Right */ + + (art[left] + art[top]) * x[top_left] /* Top Left */ + - (art[right] + art[top]) * x[top_right] /* Top Right */ + ); if (i_r == grid.numberSmootherCircles()) { const double coeff1 = 0.5 * (k1 + k2) / h1; @@ -120,7 +178,8 @@ template void SmootherTake::applyAscOrthoBlackCircleSection(ConstVector x, ConstVector rhs, Vector temp) { - using smoother_take::nodeApplyAscOrthoCircleTake; + using smoother_take::nodeApplyAscOrthoCircleTakeBoundary; + using smoother_take::nodeApplyAscOrthoCircleTakeInterior; const PolarGrid& grid = Smoother::grid_; const LevelCacheType& level_cache = Smoother::level_cache_; @@ -136,17 +195,44 @@ void SmootherTake::applyAscOrthoBlackCircleSection(ConstVector>( // Rank of the index space - {0, 0}, // Starting point of the index space - {num_black_circles, grid.ntheta()} // Ending point of the index space - ), - // Kokkos lambda function to execute for each point in the index space - KOKKOS_LAMBDA(const int circle_task, const int i_theta) { - int i_r = start_black_circles + circle_task * 2; - nodeApplyAscOrthoCircleTake(i_r, i_theta, grid, DirBC_Interior, x, rhs, temp, arr, att, art); - }); + if (start_black_circles == 0) { + Kokkos::parallel_for( + "Smoother Take: ApplyAscOrtho (Black Circular)", + Kokkos::MDRangePolicy>( // Rank of the index space + {0, 0}, // Starting point of the index space + {1, grid.ntheta()} // Ending point of the index space + ), + // Kokkos lambda function to execute for each point in the index space + KOKKOS_LAMBDA(const int circle_task, const int i_theta) { + int i_r = start_black_circles + circle_task * 2; + nodeApplyAscOrthoCircleTakeBoundary(i_r, i_theta, grid, DirBC_Interior, x, rhs, temp, arr, att, art); + }); + + Kokkos::parallel_for( + "Smoother Take: ApplyAscOrtho (Black Circular)", + Kokkos::MDRangePolicy>( // Rank of the index space + {1, 0}, // Starting point of the index space + {num_black_circles, grid.ntheta()} // Ending point of the index space + ), + // Kokkos lambda function to execute for each point in the index space + KOKKOS_LAMBDA(const int circle_task, const int i_theta) { + int i_r = start_black_circles + circle_task * 2; + nodeApplyAscOrthoCircleTakeInterior(i_r, i_theta, grid, DirBC_Interior, x, rhs, temp, arr, att, art); + }); + } + else { + Kokkos::parallel_for( + "Smoother Take: ApplyAscOrtho (Black Circular)", + Kokkos::MDRangePolicy>( // Rank of the index space + {0, 0}, // Starting point of the index space + {num_black_circles, grid.ntheta()} // Ending point of the index space + ), + // Kokkos lambda function to execute for each point in the index space + KOKKOS_LAMBDA(const int circle_task, const int i_theta) { + int i_r = start_black_circles + circle_task * 2; + nodeApplyAscOrthoCircleTakeInterior(i_r, i_theta, grid, DirBC_Interior, x, rhs, temp, arr, att, art); + }); + } Kokkos::fence(); } @@ -155,7 +241,8 @@ template void SmootherTake::applyAscOrthoWhiteCircleSection(ConstVector x, ConstVector rhs, Vector temp) { - using smoother_take::nodeApplyAscOrthoCircleTake; + using smoother_take::nodeApplyAscOrthoCircleTakeBoundary; + using smoother_take::nodeApplyAscOrthoCircleTakeInterior; const PolarGrid& grid = Smoother::grid_; const LevelCacheType& level_cache = Smoother::level_cache_; @@ -171,17 +258,44 @@ void SmootherTake::applyAscOrthoWhiteCircleSection(ConstVector>( // Rank of the index space - {0, 0}, // Starting point of the index space - {num_white_circles, grid.ntheta()} // Ending point of the index space - ), - // Kokkos lambda function to execute for each point in the index space - KOKKOS_LAMBDA(const int circle_task, const int i_theta) { - const int i_r = start_white_circles + circle_task * 2; - nodeApplyAscOrthoCircleTake(i_r, i_theta, grid, DirBC_Interior, x, rhs, temp, arr, att, art); - }); + if (start_white_circles == 0) { + Kokkos::parallel_for( + "Smoother Take: ApplyAscOrtho (White Circular)", + Kokkos::MDRangePolicy>( // Rank of the index space + {0, 0}, // Starting point of the index space + {1, grid.ntheta()} // Ending point of the index space + ), + // Kokkos lambda function to execute for each point in the index space + KOKKOS_LAMBDA(const int circle_task, const int i_theta) { + const int i_r = start_white_circles + circle_task * 2; + nodeApplyAscOrthoCircleTakeBoundary(i_r, i_theta, grid, DirBC_Interior, x, rhs, temp, arr, att, art); + }); + + Kokkos::parallel_for( + "Smoother Take: ApplyAscOrtho (White Circular)", + Kokkos::MDRangePolicy>( // Rank of the index space + {1, 0}, // Starting point of the index space + {num_white_circles, grid.ntheta()} // Ending point of the index space + ), + // Kokkos lambda function to execute for each point in the index space + KOKKOS_LAMBDA(const int circle_task, const int i_theta) { + const int i_r = start_white_circles + circle_task * 2; + nodeApplyAscOrthoCircleTakeInterior(i_r, i_theta, grid, DirBC_Interior, x, rhs, temp, arr, att, art); + }); + } + else { + Kokkos::parallel_for( + "Smoother Take: ApplyAscOrtho (White Circular)", + Kokkos::MDRangePolicy>( // Rank of the index space + {0, 0}, // Starting point of the index space + {num_white_circles, grid.ntheta()} // Ending point of the index space + ), + // Kokkos lambda function to execute for each point in the index space + KOKKOS_LAMBDA(const int circle_task, const int i_theta) { + const int i_r = start_white_circles + circle_task * 2; + nodeApplyAscOrthoCircleTakeInterior(i_r, i_theta, grid, DirBC_Interior, x, rhs, temp, arr, att, art); + }); + } Kokkos::fence(); } @@ -190,7 +304,8 @@ template void SmootherTake::applyAscOrthoBlackRadialSection(ConstVector x, ConstVector rhs, Vector temp) { - using smoother_take::nodeApplyAscOrthoRadialTake; + using smoother_take::nodeApplyAscOrthoRadialTakeBoundary; + using smoother_take::nodeApplyAscOrthoRadialTakeInterior; const PolarGrid& grid = Smoother::grid_; const LevelCacheType& level_cache = Smoother::level_cache_; @@ -210,12 +325,36 @@ void SmootherTake::applyAscOrthoBlackRadialSection(ConstVector>( // Rank of the index space {0, grid.numberSmootherCircles()}, // Starting point of the index space + {num_black_radial_lines, grid.numberSmootherCircles() + 1} // Ending point of the index space + ), + // Kokkos lambda function to execute for each point in the index space + KOKKOS_LAMBDA(const int radial_task, const int i_r) { + const int i_theta = start_black_radials + radial_task * 2; + nodeApplyAscOrthoRadialTakeBoundary(i_r, i_theta, grid, DirBC_Interior, x, rhs, temp, arr, att, art); + }); + + Kokkos::parallel_for( + "Smoother Take: ApplyAscOrtho (Black Radial)", + Kokkos::MDRangePolicy>( // Rank of the index space + {0, grid.numberSmootherCircles() + 1}, // Starting point of the index space + {num_black_radial_lines, grid.nr() - 2} // Ending point of the index space + ), + // Kokkos lambda function to execute for each point in the index space + KOKKOS_LAMBDA(const int radial_task, const int i_r) { + const int i_theta = start_black_radials + radial_task * 2; + nodeApplyAscOrthoRadialTakeInterior(i_r, i_theta, grid, DirBC_Interior, x, rhs, temp, arr, att, art); + }); + + Kokkos::parallel_for( + "Smoother Take: ApplyAscOrtho (Black Radial)", + Kokkos::MDRangePolicy>( // Rank of the index space + {0, grid.nr() - 2}, // Starting point of the index space {num_black_radial_lines, grid.nr()} // Ending point of the index space ), // Kokkos lambda function to execute for each point in the index space KOKKOS_LAMBDA(const int radial_task, const int i_r) { const int i_theta = start_black_radials + radial_task * 2; - nodeApplyAscOrthoRadialTake(i_r, i_theta, grid, DirBC_Interior, x, rhs, temp, arr, att, art); + nodeApplyAscOrthoRadialTakeBoundary(i_r, i_theta, grid, DirBC_Interior, x, rhs, temp, arr, att, art); }); Kokkos::fence(); @@ -225,7 +364,8 @@ template void SmootherTake::applyAscOrthoWhiteRadialSection(ConstVector x, ConstVector rhs, Vector temp) { - using smoother_take::nodeApplyAscOrthoRadialTake; + using smoother_take::nodeApplyAscOrthoRadialTakeBoundary; + using smoother_take::nodeApplyAscOrthoRadialTakeInterior; const PolarGrid& grid = Smoother::grid_; const LevelCacheType& level_cache = Smoother::level_cache_; @@ -245,12 +385,36 @@ void SmootherTake::applyAscOrthoWhiteRadialSection(ConstVector>( // Rank of the index space {0, grid.numberSmootherCircles()}, // Starting point of the index space + {num_white_radial_lines, grid.numberSmootherCircles() + 1} // Ending point of the index space + ), + // Kokkos lambda function to execute for each point in the index space + KOKKOS_LAMBDA(const int radial_task, const int i_r) { + const int i_theta = start_white_radials + radial_task * 2; + nodeApplyAscOrthoRadialTakeBoundary(i_r, i_theta, grid, DirBC_Interior, x, rhs, temp, arr, att, art); + }); + + Kokkos::parallel_for( + "Smoother Take: ApplyAscOrtho (White Radial)", + Kokkos::MDRangePolicy>( // Rank of the index space + {0, grid.numberSmootherCircles() + 1}, // Starting point of the index space + {num_white_radial_lines, grid.nr() - 2} // Ending point of the index space + ), + // Kokkos lambda function to execute for each point in the index space + KOKKOS_LAMBDA(const int radial_task, const int i_r) { + const int i_theta = start_white_radials + radial_task * 2; + nodeApplyAscOrthoRadialTakeInterior(i_r, i_theta, grid, DirBC_Interior, x, rhs, temp, arr, att, art); + }); + + Kokkos::parallel_for( + "Smoother Take: ApplyAscOrtho (White Radial)", + Kokkos::MDRangePolicy>( // Rank of the index space + {0, grid.nr() - 2}, // Starting point of the index space {num_white_radial_lines, grid.nr()} // Ending point of the index space ), // Kokkos lambda function to execute for each point in the index space KOKKOS_LAMBDA(const int radial_task, const int i_r) { const int i_theta = start_white_radials + radial_task * 2; - nodeApplyAscOrthoRadialTake(i_r, i_theta, grid, DirBC_Interior, x, rhs, temp, arr, att, art); + nodeApplyAscOrthoRadialTakeBoundary(i_r, i_theta, grid, DirBC_Interior, x, rhs, temp, arr, att, art); }); Kokkos::fence(); From 96865a0468cc62868a164e7057aed1275fe64bfc Mon Sep 17 00:00:00 2001 From: julianlitz Date: Mon, 17 Aug 2026 00:17:14 +0200 Subject: [PATCH 2/7] Add ExtrapolatedSmoother --- .../applyAscOrtho.inl | 184 +++++++++++++++++- 1 file changed, 182 insertions(+), 2 deletions(-) diff --git a/include/ExtrapolatedSmoother/ExtrapolatedSmootherTake/applyAscOrtho.inl b/include/ExtrapolatedSmoother/ExtrapolatedSmootherTake/applyAscOrtho.inl index bc8bc5ca..b4dbb66d 100644 --- a/include/ExtrapolatedSmoother/ExtrapolatedSmootherTake/applyAscOrtho.inl +++ b/include/ExtrapolatedSmoother/ExtrapolatedSmootherTake/applyAscOrtho.inl @@ -3,6 +3,110 @@ namespace extrapolated_smoother_take { +static KOKKOS_INLINE_FUNCTION void +nodeApplyAscOrthoCircleTakeInterior(const int i_r, const int i_theta, const PolarGrid& grid, const bool DirBC_Interior, + ConstVector& x, ConstVector& rhs, Vector& result, + ConstVector& arr, ConstVector& att, ConstVector& art) +{ + KOKKOS_ASSERT(i_r >= 1 && i_r <= grid.numberSmootherCircles()); + + if (!(i_r & 1) && !(i_theta & 1)) { + const int center = grid.index(i_r, i_theta); + result[center] = x[center]; + return; + } + + const int i_theta_M1 = grid.wrapThetaIndex(i_theta - 1); + const int i_theta_P1 = grid.wrapThetaIndex(i_theta + 1); + + const int bottom_left = grid.index(i_r - 1, i_theta_M1); + const int left = grid.index(i_r - 1, i_theta); + const int top_left = grid.index(i_r - 1, i_theta_P1); + 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); + const int bottom_right = grid.index(i_r + 1, i_theta_M1); + const int right = grid.index(i_r + 1, i_theta); + const int top_right = grid.index(i_r + 1, i_theta_P1); + + 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; + + result[center] = rhs[center] - (-coeff1 * (arr[center] + arr[left]) * x[left] /* Left */ + - coeff2 * (arr[center] + arr[right]) * x[right] /* Right */ + + - (art[left] + art[bottom]) * x[bottom_left] /* Bottom Left */ + + (art[right] + art[bottom]) * x[bottom_right] /* Bottom Right */ + + (art[left] + art[top]) * x[top_left] /* Top Left */ + - (art[right] + art[top]) * x[top_right] /* Top Right */ + ); + + if (!(i_r & 1) && (i_theta & 1)) { + result[center] -= (-coeff3 * (att[center] + att[bottom]) * x[bottom] /* Bottom */ + - coeff4 * (att[center] + att[top]) * x[top] /* Top */ + ); + } +} + +static KOKKOS_INLINE_FUNCTION void +nodeApplyAscOrthoRadialTakeInterior(const int i_r, const int i_theta, const PolarGrid& grid, const bool DirBC_Interior, + ConstVector& x, ConstVector& rhs, Vector& result, + ConstVector& arr, ConstVector& att, ConstVector& art) +{ + assert(i_r > grid.numberSmootherCircles() - 1 && i_r < grid.nr() - 2); + + if (!(i_r & 1) && !(i_theta & 1)) { + const int center = grid.index(i_r, i_theta); + result[center] = x[center]; + return; + } + + const int i_theta_M1 = grid.wrapThetaIndex(i_theta - 1); + const int i_theta_P1 = grid.wrapThetaIndex(i_theta + 1); + + const int bottom_left = grid.index(i_r - 1, i_theta_M1); + const int bottom = grid.index(i_r, i_theta_M1); + const int bottom_right = grid.index(i_r + 1, i_theta_M1); + const int left = grid.index(i_r - 1, i_theta); + const int center = grid.index(i_r, i_theta); + const int right = grid.index(i_r + 1, i_theta); + const int top_left = grid.index(i_r - 1, i_theta_P1); + const int top = grid.index(i_r, i_theta_P1); + const int top_right = grid.index(i_r + 1, i_theta_P1); + + const double h1 = grid.radialSpacing(i_r - 1); + const double h2 = grid.radialSpacing(i_r); + const double k1 = grid.angularSpacing(i_theta - 1); + const double k2 = grid.angularSpacing(i_theta); + + const double coeff1 = 0.5 * (k1 + k2) / h1; + const double coeff2 = 0.5 * (k1 + k2) / h2; + const double coeff3 = 0.5 * (h1 + h2) / k1; + const double coeff4 = 0.5 * (h1 + h2) / k2; + + result[center] = rhs[center] - (-coeff3 * (att[center] + att[bottom]) * x[bottom] /* Bottom */ + - coeff4 * (att[center] + att[top]) * x[top] /* Top */ + + - (art[left] + art[bottom]) * x[bottom_left] /* Bottom Left */ + + (art[right] + art[bottom]) * x[bottom_right] /* Bottom Right */ + + (art[left] + art[top]) * x[top_left] /* Top Left */ + - (art[right] + art[top]) * x[top_right] /* Top Right */ + ); + + if (!(i_theta & 1) && (i_r & 1)) { + result[center] -= (-coeff1 * (arr[center] + arr[left]) * x[left] /* Left */ + - coeff2 * (arr[center] + arr[right]) * x[right] /* Right */ + ); + } +} + static KOKKOS_INLINE_FUNCTION void nodeApplyAscOrthoCircleTake(const int i_r, const int i_theta, const PolarGrid& grid, const bool DirBC_Interior, ConstVector& x, ConstVector& rhs, Vector& result, @@ -465,6 +569,7 @@ void ExtrapolatedSmootherTake::applyAscOrthoBlackCircleSection(C Vector temp) { using extrapolated_smoother_take::nodeApplyAscOrthoCircleTake; + using extrapolated_smoother_take::nodeApplyAscOrthoCircleTakeInterior; const PolarGrid& grid = ExtrapolatedSmoother::grid_; const LevelCacheType& level_cache = ExtrapolatedSmoother::level_cache_; @@ -484,7 +589,7 @@ void ExtrapolatedSmootherTake::applyAscOrthoBlackCircleSection(C "ExtrapolatedSmootherTake: ApplyAscOrtho (Black Circular)", Kokkos::MDRangePolicy>( // Rank of the index space {0, 0}, // Starting point of the index space - {num_black_circles, grid.ntheta()} // Ending point of the index space + {1, grid.ntheta()} // Ending point of the index space ), // Kokkos lambda function to execute for each point in the index space KOKKOS_LAMBDA(const int circle_task, const int i_theta) { @@ -492,6 +597,18 @@ void ExtrapolatedSmootherTake::applyAscOrthoBlackCircleSection(C nodeApplyAscOrthoCircleTake(i_r, i_theta, grid, DirBC_Interior, x, rhs, temp, arr, att, art); }); + Kokkos::parallel_for( + "ExtrapolatedSmootherTake: ApplyAscOrtho (Black Circular)", + Kokkos::MDRangePolicy>( // Rank of the index space + {1, 0}, // Starting point of the index space + {num_black_circles, grid.ntheta()} // Ending point of the index space + ), + // Kokkos lambda function to execute for each point in the index space + KOKKOS_LAMBDA(const int circle_task, const int i_theta) { + int i_r = start_black_circles + circle_task * 2; + nodeApplyAscOrthoCircleTakeInterior(i_r, i_theta, grid, DirBC_Interior, x, rhs, temp, arr, att, art); + }); + Kokkos::fence(); } @@ -501,6 +618,7 @@ void ExtrapolatedSmootherTake::applyAscOrthoWhiteCircleSection(C Vector temp) { using extrapolated_smoother_take::nodeApplyAscOrthoCircleTake; + using extrapolated_smoother_take::nodeApplyAscOrthoCircleTakeInterior; const PolarGrid& grid = ExtrapolatedSmoother::grid_; const LevelCacheType& level_cache = ExtrapolatedSmoother::level_cache_; @@ -520,7 +638,7 @@ void ExtrapolatedSmootherTake::applyAscOrthoWhiteCircleSection(C "ExtrapolatedSmootherTake: ApplyAscOrtho (White Circular)", Kokkos::MDRangePolicy>( // Rank of the index space {0, 0}, // Starting point of the index space - {num_white_circles, grid.ntheta()} // Ending point of the index space + {1, grid.ntheta()} // Ending point of the index space ), // Kokkos lambda function to execute for each point in the index space KOKKOS_LAMBDA(const int circle_task, const int i_theta) { @@ -528,6 +646,18 @@ void ExtrapolatedSmootherTake::applyAscOrthoWhiteCircleSection(C nodeApplyAscOrthoCircleTake(i_r, i_theta, grid, DirBC_Interior, x, rhs, temp, arr, att, art); }); + Kokkos::parallel_for( + "ExtrapolatedSmootherTake: ApplyAscOrtho (White Circular)", + Kokkos::MDRangePolicy>( // Rank of the index space + {1, 0}, // Starting point of the index space + {num_white_circles, grid.ntheta()} // Ending point of the index space + ), + // Kokkos lambda function to execute for each point in the index space + KOKKOS_LAMBDA(const int circle_task, const int i_theta) { + const int i_r = start_white_circles + circle_task * 2; + nodeApplyAscOrthoCircleTakeInterior(i_r, i_theta, grid, DirBC_Interior, x, rhs, temp, arr, att, art); + }); + Kokkos::fence(); } @@ -537,6 +667,7 @@ void ExtrapolatedSmootherTake::applyAscOrthoBlackRadialSection(C Vector temp) { using extrapolated_smoother_take::nodeApplyAscOrthoRadialTake; + using extrapolated_smoother_take::nodeApplyAscOrthoRadialTakeInterior; const PolarGrid& grid = ExtrapolatedSmoother::grid_; const LevelCacheType& level_cache = ExtrapolatedSmoother::level_cache_; @@ -556,6 +687,30 @@ void ExtrapolatedSmootherTake::applyAscOrthoBlackRadialSection(C "ExtrapolatedSmootherTake: ApplyAscOrtho (Black Radial)", Kokkos::MDRangePolicy>( // Rank of the index space {0, grid.numberSmootherCircles()}, // Starting point of the index space + {num_black_radial_lines, grid.numberSmootherCircles() + 1} // Ending point of the index space + ), + // Kokkos lambda function to execute for each point in the index space + KOKKOS_LAMBDA(const int radial_task, const int i_r) { + const int i_theta = start_black_radials + radial_task * 2; + nodeApplyAscOrthoRadialTake(i_r, i_theta, grid, DirBC_Interior, x, rhs, temp, arr, att, art); + }); + + Kokkos::parallel_for( + "ExtrapolatedSmootherTake: ApplyAscOrtho (Black Radial)", + Kokkos::MDRangePolicy>( // Rank of the index space + {0, grid.numberSmootherCircles() + 1}, // Starting point of the index space + {num_black_radial_lines, grid.nr() - 2} // Ending point of the index space + ), + // Kokkos lambda function to execute for each point in the index space + KOKKOS_LAMBDA(const int radial_task, const int i_r) { + const int i_theta = start_black_radials + radial_task * 2; + nodeApplyAscOrthoRadialTakeInterior(i_r, i_theta, grid, DirBC_Interior, x, rhs, temp, arr, att, art); + }); + + Kokkos::parallel_for( + "ExtrapolatedSmootherTake: ApplyAscOrtho (Black Radial)", + Kokkos::MDRangePolicy>( // Rank of the index space + {0, grid.nr() - 2}, // Starting point of the index space {num_black_radial_lines, grid.nr()} // Ending point of the index space ), // Kokkos lambda function to execute for each point in the index space @@ -573,6 +728,7 @@ void ExtrapolatedSmootherTake::applyAscOrthoWhiteRadialSection(C Vector temp) { using extrapolated_smoother_take::nodeApplyAscOrthoRadialTake; + using extrapolated_smoother_take::nodeApplyAscOrthoRadialTakeInterior; const PolarGrid& grid = ExtrapolatedSmoother::grid_; const LevelCacheType& level_cache = ExtrapolatedSmoother::level_cache_; @@ -592,6 +748,30 @@ void ExtrapolatedSmootherTake::applyAscOrthoWhiteRadialSection(C "ExtrapolatedSmootherTake: ApplyAscOrtho (White Radial)", Kokkos::MDRangePolicy>( // Rank of the index space {0, grid.numberSmootherCircles()}, // Starting point of the index space + {num_white_radial_lines, grid.numberSmootherCircles() + 1} // Ending point of the index space + ), + // Kokkos lambda function to execute for each point in the index space + KOKKOS_LAMBDA(const int radial_task, const int i_r) { + const int i_theta = start_white_radials + radial_task * 2; + nodeApplyAscOrthoRadialTake(i_r, i_theta, grid, DirBC_Interior, x, rhs, temp, arr, att, art); + }); + + Kokkos::parallel_for( + "ExtrapolatedSmootherTake: ApplyAscOrtho (White Radial)", + Kokkos::MDRangePolicy>( // Rank of the index space + {0, grid.numberSmootherCircles() + 1}, // Starting point of the index space + {num_white_radial_lines, grid.nr() - 2} // Ending point of the index space + ), + // Kokkos lambda function to execute for each point in the index space + KOKKOS_LAMBDA(const int radial_task, const int i_r) { + const int i_theta = start_white_radials + radial_task * 2; + nodeApplyAscOrthoRadialTakeInterior(i_r, i_theta, grid, DirBC_Interior, x, rhs, temp, arr, att, art); + }); + + Kokkos::parallel_for( + "ExtrapolatedSmootherTake: ApplyAscOrtho (White Radial)", + Kokkos::MDRangePolicy>( // Rank of the index space + {0, grid.nr() - 2}, // Starting point of the index space {num_white_radial_lines, grid.nr()} // Ending point of the index space ), // Kokkos lambda function to execute for each point in the index space From c875ba78a6147c721fe33e721447f4ac9bda1580 Mon Sep 17 00:00:00 2001 From: julianlitz Date: Mon, 17 Aug 2026 00:22:50 +0200 Subject: [PATCH 3/7] Format --- include/Residual/ResidualTake/applyATake.inl | 32 ++++++++++---------- 1 file changed, 16 insertions(+), 16 deletions(-) diff --git a/include/Residual/ResidualTake/applyATake.inl b/include/Residual/ResidualTake/applyATake.inl index daeb5b12..bef109e4 100644 --- a/include/Residual/ResidualTake/applyATake.inl +++ b/include/Residual/ResidualTake/applyATake.inl @@ -36,17 +36,17 @@ static KOKKOS_INLINE_FUNCTION void applySystemMatrixTakeInterior(const int i_r, const double coeff4 = 0.5 * (h1 + h2) / k2; const double coeff5 = 0.25 * (h1 + h2) * (k1 + k2); - result[center] = (+coeff5 * coeff_beta[center] * Kokkos::fabs(detDF[center]) * x[center]; /* beta_{i,j} */ + result[center] = (+coeff5 * coeff_beta[center] * Kokkos::fabs(detDF[center]) * x[center] /* beta_{i,j} */ - +coeff1 * (arr[center] + arr[left]) * (x[center] - x[left]); /* Center: (Left) - Left */ - +coeff2 * (arr[center] + arr[right]) * (x[center] - x[right]); /* Center: (Right) - Right */ - +coeff3 * (att[center] + att[bottom]) * (x[center] - x[bottom]); /* Center: (Bottom) - Bottom */ - +coeff4 * (att[center] + att[top]) * (x[center] - x[top]); /* Center: (Top) - Top */ + + coeff1 * (arr[center] + arr[left]) * (x[center] - x[left]) /* Center: (Left) - Left */ + + coeff2 * (arr[center] + arr[right]) * (x[center] - x[right]) /* Center: (Right) - Right */ + + coeff3 * (att[center] + att[bottom]) * (x[center] - x[bottom]) /* Center: (Bottom) - Bottom */ + + coeff4 * (att[center] + att[top]) * (x[center] - x[top]) /* Center: (Top) - Top */ - -(art[left] + art[bottom]) * x[bottom_left]; /* Bottom Left */ - +(art[left] + art[top]) * x[top_left]; /* Top Left */ - +(art[right] + art[bottom]) * x[bottom_right]; /* Bottom Right */ - -(art[right] + art[top]) * x[top_right]; /* Top Right */ + - (art[left] + art[bottom]) * x[bottom_left] /* Bottom Left */ + + (art[left] + art[top]) * x[top_left] /* Top Left */ + + (art[right] + art[bottom]) * x[bottom_right] /* Bottom Right */ + - (art[right] + art[top]) * x[top_right] /* Top Right */ ); } @@ -87,15 +87,15 @@ static KOKKOS_INLINE_FUNCTION void applySystemMatrixTakeBoundary(const int i_r, const double coeff4 = 0.5 * (h1 + h2) / k2; const double coeff5 = 0.25 * (h1 + h2) * (k1 + k2); - result[center] = (+coeff5 * coeff_beta[center] * Kokkos::fabs(detDF[center]) * x[center]; /* beta_{i,j} */ + result[center] = (+coeff5 * coeff_beta[center] * Kokkos::fabs(detDF[center]) * x[center] /* beta_{i,j} */ - +coeff1 * (arr[center] + arr[left]) * (x[center] - x[left]); /* Center: (Left) - Left */ - +coeff2 * (arr[center] + arr[right]) * (x[center] - x[right]); /* Center: (Right) - Right */ - +coeff3 * (att[center] + att[bottom]) * (x[center] - x[bottom]); /* Center: (Bottom) - Bottom */ - +coeff4 * (att[center] + att[top]) * (x[center] - x[top]); /* Center: (Top) - Top */ + + coeff1 * (arr[center] + arr[left]) * (x[center] - x[left]) /* Center: (Left) - Left */ + + coeff2 * (arr[center] + arr[right]) * (x[center] - x[right]) /* Center: (Right) - Right */ + + coeff3 * (att[center] + att[bottom]) * (x[center] - x[bottom]) /* Center: (Bottom) - Bottom */ + + coeff4 * (att[center] + att[top]) * (x[center] - x[top]) /* Center: (Top) - Top */ - +(art[right] + art[bottom]) * x[bottom_right]; /* Bottom Right */ - -(art[right] + art[top]) * x[top_right]; /* Top Right */ + + (art[right] + art[bottom]) * x[bottom_right] /* Bottom Right */ + - (art[right] + art[top]) * x[top_right] /* Top Right */ ); } From 4e8c5455c6f22223c4a31820453fac8becd8031a Mon Sep 17 00:00:00 2001 From: julianlitz Date: Mon, 17 Aug 2026 00:29:16 +0200 Subject: [PATCH 4/7] Fix --- include/Definitions/geometry_helper.h | 2 +- include/Residual/ResidualTake/applyATake.inl | 12 +++++++----- 2 files changed, 8 insertions(+), 6 deletions(-) diff --git a/include/Definitions/geometry_helper.h b/include/Definitions/geometry_helper.h index 5f04987b..67aa739b 100644 --- a/include/Definitions/geometry_helper.h +++ b/include/Definitions/geometry_helper.h @@ -33,7 +33,7 @@ KOKKOS_INLINE_FUNCTION void compute_jacobian_elements(const DomainGeometry& doma /* 1.0 / det(DF) * */ /* [Jtt, -Jrt] */ /* [-Jtr, Jrr] */ - /* The stencil coefficient art is scaled by 0.25 to account for the + /* The stencil coefficient art is scaled by 0.25 to account for the */ /* 9-point stencil used in the finite difference approximation */ art *= 0.25; } diff --git a/include/Residual/ResidualTake/applyATake.inl b/include/Residual/ResidualTake/applyATake.inl index bef109e4..d6f589fc 100644 --- a/include/Residual/ResidualTake/applyATake.inl +++ b/include/Residual/ResidualTake/applyATake.inl @@ -69,11 +69,13 @@ static KOKKOS_INLINE_FUNCTION void applySystemMatrixTakeBoundary(const int i_r, 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 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); - const int right = grid.index(i_r + 1, i_theta); + 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); + const int bottom_right = grid.index(i_r + 1, i_theta_M1); + const int right = grid.index(i_r + 1, i_theta); + const int top_right = grid.index(i_r + 1, i_theta_P1); // Across origin: h1 gets replaced with 2 * R0. const double h1 = 2.0 * grid.radius(0); From 75ed92f4949414e1da45534d09b8243de029f52a Mon Sep 17 00:00:00 2001 From: julianlitz Date: Sat, 22 Aug 2026 21:40:26 +0200 Subject: [PATCH 5/7] Imporve --- .../applyAscOrtho.inl | 497 +++++------------- include/Residual/ResidualGive/residualGive.h | 1 - include/Residual/ResidualTake/applyATake.inl | 18 +- include/Residual/ResidualTake/residualTake.h | 2 - .../Smoother/SmootherTake/applyAscOrtho.inl | 42 +- include/Smoother/SmootherTake/smootherTake.h | 2 - 6 files changed, 160 insertions(+), 402 deletions(-) diff --git a/include/ExtrapolatedSmoother/ExtrapolatedSmootherTake/applyAscOrtho.inl b/include/ExtrapolatedSmoother/ExtrapolatedSmootherTake/applyAscOrtho.inl index b4dbb66d..1481bbfc 100644 --- a/include/ExtrapolatedSmoother/ExtrapolatedSmootherTake/applyAscOrtho.inl +++ b/include/ExtrapolatedSmoother/ExtrapolatedSmootherTake/applyAscOrtho.inl @@ -4,11 +4,11 @@ namespace extrapolated_smoother_take { static KOKKOS_INLINE_FUNCTION void -nodeApplyAscOrthoCircleTakeInterior(const int i_r, const int i_theta, const PolarGrid& grid, const bool DirBC_Interior, - ConstVector& x, ConstVector& rhs, Vector& result, - ConstVector& arr, ConstVector& att, ConstVector& art) +nodeApplyAscOrthoCircleTakeInterior(const int i_r, const int i_theta, const PolarGrid& grid, ConstVector& x, + ConstVector& rhs, Vector& result, ConstVector& arr, + ConstVector& att, ConstVector& art) { - KOKKOS_ASSERT(i_r >= 1 && i_r <= grid.numberSmootherCircles()); + KOKKOS_ASSERT(0 < i_r && i_r <= grid.numberSmootherCircles()); if (!(i_r & 1) && !(i_theta & 1)) { const int center = grid.index(i_r, i_theta); @@ -49,6 +49,11 @@ nodeApplyAscOrthoCircleTakeInterior(const int i_r, const int i_theta, const Pola ); if (!(i_r & 1) && (i_theta & 1)) { + /* | o | x | o | */ + /* | | | | */ + /* | o | O | o | */ + /* | | | | */ + /* | o | x | o | */ result[center] -= (-coeff3 * (att[center] + att[bottom]) * x[bottom] /* Bottom */ - coeff4 * (att[center] + att[top]) * x[top] /* Top */ ); @@ -56,11 +61,58 @@ nodeApplyAscOrthoCircleTakeInterior(const int i_r, const int i_theta, const Pola } static KOKKOS_INLINE_FUNCTION void -nodeApplyAscOrthoRadialTakeInterior(const int i_r, const int i_theta, const PolarGrid& grid, const bool DirBC_Interior, +nodeApplyAscOrthoCircleTakeBoundary(const int i_r, const int i_theta, const PolarGrid& grid, const bool DirBC_Interior, ConstVector& x, ConstVector& rhs, Vector& result, ConstVector& arr, ConstVector& att, ConstVector& art) { - assert(i_r > grid.numberSmootherCircles() - 1 && i_r < grid.nr() - 2); + KOKKOS_ASSERT(i_r == 0); + + if (!(i_theta & 1)) { + const int center = grid.index(i_r, i_theta); + result[center] = x[center]; + return; + } + + if (DirBC_Interior) { + const int center = grid.index(i_r, i_theta); + result[center] = rhs[center]; + } + else { + const double h1 = 2.0 * grid.radius(0); + const double h2 = grid.radialSpacing(i_r); + const double k1 = grid.angularSpacing(i_theta - 1); + const double k2 = grid.angularSpacing(i_theta); + + const double coeff2 = 0.5 * (k1 + k2) / h2; + const 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 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); + const int bottom_right = grid.index(i_r + 1, i_theta_M1); + const int right = grid.index(i_r + 1, i_theta); + const int top_right = grid.index(i_r + 1, i_theta_P1); + + result[center] = rhs[center] - (-coeff2 * (arr[center] + arr[right]) * x[right] /* Right */ + - coeff3 * (att[center] + att[bottom]) * x[bottom] /* Bottom */ + - coeff4 * (att[center] + att[top]) * x[top] /* Top */ + + + (art[right] + art[bottom]) * x[bottom_right] /* Bottom Right */ + - (art[right] + art[top]) * x[top_right] /* Top Right */ + ); + } +} + +static KOKKOS_INLINE_FUNCTION void +nodeApplyAscOrthoRadialTakeInterior(const int i_r, const int i_theta, const PolarGrid& grid, ConstVector& x, + ConstVector& rhs, Vector& result, ConstVector& arr, + ConstVector& att, ConstVector& art) +{ + KOKKOS_ASSERT(grid.numberSmootherCircles() < i_r && i_r < grid.nr() - 2); if (!(i_r & 1) && !(i_theta & 1)) { const int center = grid.index(i_r, i_theta); @@ -101,368 +153,83 @@ nodeApplyAscOrthoRadialTakeInterior(const int i_r, const int i_theta, const Pola ); if (!(i_theta & 1) && (i_r & 1)) { + /* ---------- */ + /* o o o */ + /* ---------- */ + /* x O x */ + /* ---------- */ + /* o o o */ + /* ---------- */ result[center] -= (-coeff1 * (arr[center] + arr[left]) * x[left] /* Left */ - coeff2 * (arr[center] + arr[right]) * x[right] /* Right */ ); } } -static KOKKOS_INLINE_FUNCTION void nodeApplyAscOrthoCircleTake(const int i_r, const int i_theta, const PolarGrid& grid, - const bool DirBC_Interior, ConstVector& x, - ConstVector& rhs, Vector& result, - ConstVector& arr, ConstVector& att, - ConstVector& art) +static KOKKOS_INLINE_FUNCTION void +nodeApplyAscOrthoRadialTakeInnerBoundary(const int i_r, const int i_theta, const PolarGrid& grid, + ConstVector& x, ConstVector& rhs, Vector& result, + ConstVector& arr, ConstVector& att, ConstVector& art) { - KOKKOS_ASSERT(i_r >= 0 && i_r <= grid.numberSmootherCircles()); - - /* -------------------- */ - /* Node in the interior */ - /* -------------------- */ - if (i_r > 0 && i_r < grid.numberSmootherCircles()) { - /* -------------------------- */ - /* Cyclic Tridiagonal Section */ - /* i_r % 2 == 1 */ - const double h1 = grid.radialSpacing(i_r - 1); - const double h2 = grid.radialSpacing(i_r); - const double k1 = grid.angularSpacing(i_theta - 1); - const double k2 = grid.angularSpacing(i_theta); - - const double coeff1 = 0.5 * (k1 + k2) / h1; - const double coeff2 = 0.5 * (k1 + k2) / h2; - const double coeff3 = 0.5 * (h1 + h2) / k1; - const double coeff4 = 0.5 * (h1 + h2) / k2; + KOKKOS_ASSERT(i_r == grid.numberSmootherCircles()); - const int i_theta_M1 = grid.wrapThetaIndex(i_theta - 1); - const int i_theta_P1 = grid.wrapThetaIndex(i_theta + 1); - - const int bottom_left = grid.index(i_r - 1, i_theta_M1); - const int left = grid.index(i_r - 1, i_theta); - const int top_left = grid.index(i_r - 1, i_theta_P1); - 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); - const int bottom_right = grid.index(i_r + 1, i_theta_M1); - const int right = grid.index(i_r + 1, i_theta); - const int top_right = grid.index(i_r + 1, i_theta_P1); - - if (i_r & 1) { - /* i_r % 2 == 1 and i_theta % 2 == 1 */ - /* | x | o | x | */ - /* | | | | */ - /* | o | O | o | */ - /* | | | | */ - /* | x | o | x | */ - /* or */ - /* i_r % 2 == 1 and i_theta % 2 == 0 */ - /* | o | o | o | */ - /* | | | | */ - /* | x | O | x | */ - /* | | | | */ - /* | o | o | o | */ - result[center] = rhs[center] - (-coeff1 * (arr[center] + arr[left]) * x[left] /* Left */ - - coeff2 * (arr[center] + arr[right]) * x[right] /* Right */ - - - (art[left] + art[bottom]) * x[bottom_left] /* Bottom Left */ - + (art[right] + art[bottom]) * x[bottom_right] /* Bottom Right */ - + (art[left] + art[top]) * x[top_left] /* Top Left */ - - (art[right] + art[top]) * x[top_right] /* Top Right */ - ); - } - else { - if (i_theta & 1) { - /* i_r % 2 == 0 and i_theta % 2 == 1 */ - /* | o | x | o | */ - /* | | | | */ - /* | o | O | o | */ - /* | | | | */ - /* | o | x | o | */ - /* Fill result(i,j) */ - result[center] = rhs[center] - (-coeff1 * (arr[center] + arr[left]) * x[left] /* Left */ - - coeff2 * (arr[center] + arr[right]) * x[right] /* Right */ - - coeff3 * (att[center] + att[bottom]) * x[bottom] /* Bottom */ - - coeff4 * (att[center] + att[top]) * x[top] /* Top */ - - - (art[left] + art[bottom]) * x[bottom_left] /* Bottom Left */ - + (art[right] + art[bottom]) * x[bottom_right] /* Bottom Right */ - + (art[left] + art[top]) * x[top_left] /* Top Left */ - - (art[right] + art[top]) * x[top_right] /* Top Right */ - ); - } - else { - /* i_r % 2 == 0 and i_theta % 2 == 0 */ - /* | o | o | o | */ - /* | | | | */ - /* | o | X | o | */ - /* | | | | */ - /* | o | o | o | */ - result[center] = x[center]; - } - } - } - /* -------------------- */ - /* Node on the boundary */ - /* -------------------- */ - else if (i_r == 0) { - /* ------------------------------------------------ */ - /* Case 1: Dirichlet boundary on the inner boundary */ - /* ------------------------------------------------ */ + if (!(i_r & 1) && !(i_theta & 1)) { const int center = grid.index(i_r, i_theta); - if (DirBC_Interior) { - if (i_theta & 1) { - /* i_theta % 2 == 1 */ - /* || x | o | x | */ - /* || | | | */ - /* || O | o | o | */ - /* || | | | */ - /* || x | o | x | */ - result[center] = rhs[center]; - } - else { - /* i_theta % 2 == 0 */ - /* || o | o | o | */ - /* || | | | */ - /* || X | o | x | */ - /* || | | | */ - /* || o | o | o | */ - result[center] = x[center]; - } - } - 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 double h1 = 2.0 * grid.radius(0); - const double h2 = grid.radialSpacing(i_r); - const double k1 = grid.angularSpacing(i_theta - 1); - const double k2 = grid.angularSpacing(i_theta); - - const double coeff2 = 0.5 * (k1 + k2) / h2; - const 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 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); - const int bottom_right = grid.index(i_r + 1, i_theta_M1); - const int right = grid.index(i_r + 1, i_theta); - const int top_right = grid.index(i_r + 1, i_theta_P1); - - if (i_theta & 1) { - /* i_theta % 2 == 1 */ - /* -| x | o | x | */ - /* -| | | | */ - /* -| O | o | o | */ - /* -| | | | */ - /* -| x | o | x | */ - result[center] = - rhs[center] - - (-coeff2 * (arr[center] + arr[right]) * x[right] /* Right */ - - coeff3 * (att[center] + att[bottom]) * x[bottom] /* Bottom */ - - coeff4 * (att[center] + att[top]) * x[top] /* Top */ - - /* - (art[left] + art[bottom]) * x[bottom_left] // Bottom Left: REMOVED DUE TO ARTIFICAL 7 POINT STENCIL */ - + (art[right] + art[bottom]) * x[bottom_right] /* Bottom Right */ - - /* + (art[left] + art[top]) * x[top_left] // Top Left: REMOVED DUE TO ARTIFICAL 7 POINT STENCIL */ - - (art[right] + art[top]) * x[top_right] /* Top Right */ - ); - } - else { - /* i_theta % 2 == 0 */ - /* -| o | o | o | */ - /* -| | | | */ - /* -| X | o | x | */ - /* -| | | | */ - /* -| o | o | o | */ - result[center] = x[center]; - } - } + result[center] = x[center]; + return; } -} -static KOKKOS_INLINE_FUNCTION void nodeApplyAscOrthoRadialTake(const int i_r, const int i_theta, const PolarGrid& grid, - const bool DirBC_Interior, ConstVector& x, - ConstVector& rhs, Vector& result, - ConstVector& arr, ConstVector& att, - ConstVector& art) -{ - assert(i_r >= grid.numberSmootherCircles() - 1 && i_r < grid.nr()); + 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); - /* -------------------- */ - /* Node in the interior */ - /* -------------------- */ - if (i_r > grid.numberSmootherCircles() && i_r < grid.nr() - 2) { - const double h1 = grid.radialSpacing(i_r - 1); - const double h2 = grid.radialSpacing(i_r); - const double k1 = grid.angularSpacing(i_theta - 1); - const double k2 = grid.angularSpacing(i_theta); + const double coeff1 = 0.5 * (k1 + k2) / h1; + const double coeff2 = 0.5 * (k1 + k2) / h2; + const double coeff3 = 0.5 * (h1 + h2) / k1; + const double coeff4 = 0.5 * (h1 + h2) / k2; - const 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_M1 = grid.wrapThetaIndex(i_theta - 1); - const int i_theta_P1 = grid.wrapThetaIndex(i_theta + 1); + const int bottom_left = grid.index(i_r - 1, i_theta_M1); + const int left = grid.index(i_r - 1, i_theta); + const int top_left = grid.index(i_r - 1, i_theta_P1); + 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); + const int bottom_right = grid.index(i_r + 1, i_theta_M1); + const int right = grid.index(i_r + 1, i_theta); + const int top_right = grid.index(i_r + 1, i_theta_P1); - const int bottom_left = grid.index(i_r - 1, i_theta_M1); - const int left = grid.index(i_r - 1, i_theta); - const int top_left = grid.index(i_r - 1, i_theta_P1); - 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); - const int bottom_right = grid.index(i_r + 1, i_theta_M1); - const int right = grid.index(i_r + 1, i_theta); - const int top_right = grid.index(i_r + 1, i_theta_P1); + result[center] = rhs[center] - (-coeff1 * (arr[center] + arr[left]) * x[left] /* Left */ + - coeff3 * (att[center] + att[bottom]) * x[bottom] /* Bottom */ + - coeff4 * (att[center] + att[top]) * x[top] /* Top */ - if (i_theta & 1) { - /* i_theta % 2 == 1 and i_r % 2 == 1 */ - /* ---------- */ - /* x o x */ - /* ---------- */ - /* o O o */ - /* ---------- */ - /* x o x */ - /* ---------- */ - /* or */ - /* i_theta % 2 == 1 and i_r % 2 == 0 */ - /* ---------- */ - /* o x o */ - /* ---------- */ - /* o O o */ - /* ---------- */ - /* o x o */ - /* ---------- */ - result[center] = rhs[center] - (-coeff3 * (att[center] + att[bottom]) * x[bottom] /* Bottom */ - - coeff4 * (att[center] + att[top]) * x[top] /* Top */ + - (art[left] + art[bottom]) * x[bottom_left] /* Bottom Left */ + + (art[right] + art[bottom]) * x[bottom_right] /* Bottom Right */ + + (art[left] + art[top]) * x[top_left] /* Top Left */ + - (art[right] + art[top]) * x[top_right] /* Top Right */ + ); - - (art[left] + art[bottom]) * x[bottom_left] /* Bottom Left */ - + (art[right] + art[bottom]) * x[bottom_right] /* Bottom Right */ - + (art[left] + art[top]) * x[top_left] /* Top Left */ - - (art[right] + art[top]) * x[top_right] /* Top Right */ - ); - } - else { - if (i_r & 1) { - /* i_theta % 2 == 0 and i_r % 2 == 1 */ - /* ---------- */ - /* o o o */ - /* ---------- */ - /* x O x */ - /* ---------- */ - /* o o o */ - /* ---------- */ - result[center] = rhs[center] - (-coeff1 * (arr[center] + arr[left]) * x[left] /* Left */ - - coeff2 * (arr[center] + arr[right]) * x[right] /* Right */ - - coeff3 * (att[center] + att[bottom]) * x[bottom] /* Bottom */ - - coeff4 * (att[center] + att[top]) * x[top] /* Top */ - - - (art[left] + art[bottom]) * x[bottom_left] /* Bottom Left */ - + (art[right] + art[bottom]) * x[bottom_right] /* Bottom Right */ - + (art[left] + art[top]) * x[top_left] /* Top Left */ - - (art[right] + art[top]) * x[top_right] /* Top Right */ - ); - } - else { - /* i_theta % 2 == 0 and i_r % 2 == 0 */ - /* ---------- */ - /* o o o */ - /* ---------- */ - /* o X o */ - /* ---------- */ - /* o o o */ - /* ---------- */ - result[center] = x[center]; - } - } + if (!(i_theta & 1) && (i_r & 1)) { + /* | o | o | o || o o o o */ + /* | | | || -------------- */ + /* | x | o | x || O x o x */ + /* | | | || -------------- */ + /* | o | o | o || o o o o */ + result[center] -= (-coeff2 * (arr[center] + arr[right]) * x[right] /* Right */); } - else if (i_r == grid.numberSmootherCircles()) { - const double h1 = grid.radialSpacing(i_r - 1); - const double h2 = grid.radialSpacing(i_r); - const double k1 = grid.angularSpacing(i_theta - 1); - const double k2 = grid.angularSpacing(i_theta); - - const double coeff1 = 0.5 * (k1 + k2) / h1; - const double coeff2 = 0.5 * (k1 + k2) / h2; - const double coeff3 = 0.5 * (h1 + h2) / k1; - const double coeff4 = 0.5 * (h1 + h2) / k2; - - const int i_theta_M1 = grid.wrapThetaIndex(i_theta - 1); - const int i_theta_P1 = grid.wrapThetaIndex(i_theta + 1); - - const int bottom_left = grid.index(i_r - 1, i_theta_M1); - const int left = grid.index(i_r - 1, i_theta); - const int top_left = grid.index(i_r - 1, i_theta_P1); - 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); - const int bottom_right = grid.index(i_r + 1, i_theta_M1); - const int right = grid.index(i_r + 1, i_theta); - const int top_right = grid.index(i_r + 1, i_theta_P1); +} - if (i_theta & 1) { - /* i_theta % 2 == 1 and i_r % 2 == 1 */ - /* | x | o | x || o x o x */ - /* | | | || -------------- */ - /* | o | o | o || O o o o */ - /* | | | || -------------- */ - /* | x | o | x || o x o x */ - /* or */ - /* i_theta % 2 == 1 and i_r % 2 == 0 */ - /* | o | x | o || x o x o */ - /* | | | || -------------- */ - /* | o | o | o || O o o o */ - /* | | | || -------------- */ - /* | o | x | o || x o x o */ - result[center] = rhs[center] - (-coeff1 * (arr[center] + arr[left]) * x[left] /* Left */ - - coeff3 * (att[center] + att[bottom]) * x[bottom] /* Bottom */ - - coeff4 * (att[center] + att[top]) * x[top] /* Top */ +static KOKKOS_INLINE_FUNCTION void +nodeApplyAscOrthoRadialTakeOuterBoundary(const int i_r, const int i_theta, const PolarGrid& grid, + ConstVector& x, ConstVector& rhs, Vector& result, + ConstVector& arr, ConstVector& att, ConstVector& art) +{ + KOKKOS_ASSERT(i_r == grid.nr() - 2 || i_r == grid.nr() - 1); - - (art[left] + art[bottom]) * x[bottom_left] /* Bottom Left */ - + (art[right] + art[bottom]) * x[bottom_right] /* Bottom Right */ - + (art[left] + art[top]) * x[top_left] /* Top Left */ - - (art[right] + art[top]) * x[top_right] /* Top Right */ - ); - } - else { - if (i_r & 1) { - /* i_theta % 2 == 0 and i_r % 2 == 1 */ - /* | o | o | o || o o o o */ - /* | | | || -------------- */ - /* | x | o | x || O x o x */ - /* | | | || -------------- */ - /* | o | o | o || o o o o */ - result[center] = rhs[center] - (-coeff1 * (arr[center] + arr[left]) * x[left] /* Left */ - - coeff2 * (arr[center] + arr[right]) * x[right] /* Right */ - - coeff3 * (att[center] + att[bottom]) * x[bottom] /* Bottom */ - - coeff4 * (att[center] + att[top]) * x[top] /* Top */ - - - (art[left] + art[bottom]) * x[bottom_left] /* Bottom Left */ - + (art[right] + art[bottom]) * x[bottom_right] /* Bottom Right */ - + (art[left] + art[top]) * x[top_left] /* Top Left */ - - (art[right] + art[top]) * x[top_right] /* Top Right */ - ); - } - else { - /* i_theta % 2 == 0 and i_r % 2 == 0 */ - /* | o | o | o || o o o o */ - /* | | | || -------------- */ - /* | o | x | o || X o x o */ - /* | | | || -------------- */ - /* | o | o | o || o o o o */ - result[center] = x[center]; - } - } - } - else if (i_r == grid.nr() - 2) { + if (i_r == grid.nr() - 2) { KOKKOS_ASSERT(i_r & 1); const double h1 = grid.radialSpacing(i_r - 1); @@ -489,7 +256,6 @@ static KOKKOS_INLINE_FUNCTION void nodeApplyAscOrthoRadialTake(const int i_r, co const int top_right = grid.index(i_r + 1, i_theta_P1); if (i_theta & 1) { - /* i_theta % 2 == 1 */ /* ---------------|| */ /* o x o x || */ /* ---------------|| */ @@ -538,7 +304,6 @@ static KOKKOS_INLINE_FUNCTION void nodeApplyAscOrthoRadialTake(const int i_r, co const int center = grid.index(i_r, i_theta); if (i_theta & 1) { - /* i_theta % 2 == 1 */ /* -----------|| */ /* x o x || */ /* -----------|| */ @@ -568,7 +333,7 @@ void ExtrapolatedSmootherTake::applyAscOrthoBlackCircleSection(C ConstVector rhs, Vector temp) { - using extrapolated_smoother_take::nodeApplyAscOrthoCircleTake; + using extrapolated_smoother_take::nodeApplyAscOrthoCircleTakeBoundary; using extrapolated_smoother_take::nodeApplyAscOrthoCircleTakeInterior; const PolarGrid& grid = ExtrapolatedSmoother::grid_; @@ -594,7 +359,7 @@ void ExtrapolatedSmootherTake::applyAscOrthoBlackCircleSection(C // Kokkos lambda function to execute for each point in the index space KOKKOS_LAMBDA(const int circle_task, const int i_theta) { int i_r = start_black_circles + circle_task * 2; - nodeApplyAscOrthoCircleTake(i_r, i_theta, grid, DirBC_Interior, x, rhs, temp, arr, att, art); + nodeApplyAscOrthoCircleTakeBoundary(i_r, i_theta, grid, DirBC_Interior, x, rhs, temp, arr, att, art); }); Kokkos::parallel_for( @@ -606,7 +371,7 @@ void ExtrapolatedSmootherTake::applyAscOrthoBlackCircleSection(C // Kokkos lambda function to execute for each point in the index space KOKKOS_LAMBDA(const int circle_task, const int i_theta) { int i_r = start_black_circles + circle_task * 2; - nodeApplyAscOrthoCircleTakeInterior(i_r, i_theta, grid, DirBC_Interior, x, rhs, temp, arr, att, art); + nodeApplyAscOrthoCircleTakeInterior(i_r, i_theta, grid, x, rhs, temp, arr, att, art); }); Kokkos::fence(); @@ -617,7 +382,7 @@ void ExtrapolatedSmootherTake::applyAscOrthoWhiteCircleSection(C ConstVector rhs, Vector temp) { - using extrapolated_smoother_take::nodeApplyAscOrthoCircleTake; + using extrapolated_smoother_take::nodeApplyAscOrthoCircleTakeBoundary; using extrapolated_smoother_take::nodeApplyAscOrthoCircleTakeInterior; const PolarGrid& grid = ExtrapolatedSmoother::grid_; @@ -643,7 +408,7 @@ void ExtrapolatedSmootherTake::applyAscOrthoWhiteCircleSection(C // Kokkos lambda function to execute for each point in the index space KOKKOS_LAMBDA(const int circle_task, const int i_theta) { const int i_r = start_white_circles + circle_task * 2; - nodeApplyAscOrthoCircleTake(i_r, i_theta, grid, DirBC_Interior, x, rhs, temp, arr, att, art); + nodeApplyAscOrthoCircleTakeBoundary(i_r, i_theta, grid, DirBC_Interior, x, rhs, temp, arr, att, art); }); Kokkos::parallel_for( @@ -655,7 +420,7 @@ void ExtrapolatedSmootherTake::applyAscOrthoWhiteCircleSection(C // Kokkos lambda function to execute for each point in the index space KOKKOS_LAMBDA(const int circle_task, const int i_theta) { const int i_r = start_white_circles + circle_task * 2; - nodeApplyAscOrthoCircleTakeInterior(i_r, i_theta, grid, DirBC_Interior, x, rhs, temp, arr, att, art); + nodeApplyAscOrthoCircleTakeInterior(i_r, i_theta, grid, x, rhs, temp, arr, att, art); }); Kokkos::fence(); @@ -666,12 +431,12 @@ void ExtrapolatedSmootherTake::applyAscOrthoBlackRadialSection(C ConstVector rhs, Vector temp) { - using extrapolated_smoother_take::nodeApplyAscOrthoRadialTake; + using extrapolated_smoother_take::nodeApplyAscOrthoRadialTakeInnerBoundary; using extrapolated_smoother_take::nodeApplyAscOrthoRadialTakeInterior; + using extrapolated_smoother_take::nodeApplyAscOrthoRadialTakeOuterBoundary; const PolarGrid& grid = ExtrapolatedSmoother::grid_; const LevelCacheType& level_cache = ExtrapolatedSmoother::level_cache_; - const bool DirBC_Interior = ExtrapolatedSmoother::DirBC_Interior_; assert(level_cache.cacheDomainGeometry()); @@ -692,7 +457,7 @@ void ExtrapolatedSmootherTake::applyAscOrthoBlackRadialSection(C // Kokkos lambda function to execute for each point in the index space KOKKOS_LAMBDA(const int radial_task, const int i_r) { const int i_theta = start_black_radials + radial_task * 2; - nodeApplyAscOrthoRadialTake(i_r, i_theta, grid, DirBC_Interior, x, rhs, temp, arr, att, art); + nodeApplyAscOrthoRadialTakeInnerBoundary(i_r, i_theta, grid, x, rhs, temp, arr, att, art); }); Kokkos::parallel_for( @@ -704,7 +469,7 @@ void ExtrapolatedSmootherTake::applyAscOrthoBlackRadialSection(C // Kokkos lambda function to execute for each point in the index space KOKKOS_LAMBDA(const int radial_task, const int i_r) { const int i_theta = start_black_radials + radial_task * 2; - nodeApplyAscOrthoRadialTakeInterior(i_r, i_theta, grid, DirBC_Interior, x, rhs, temp, arr, att, art); + nodeApplyAscOrthoRadialTakeInterior(i_r, i_theta, grid, x, rhs, temp, arr, att, art); }); Kokkos::parallel_for( @@ -716,7 +481,7 @@ void ExtrapolatedSmootherTake::applyAscOrthoBlackRadialSection(C // Kokkos lambda function to execute for each point in the index space KOKKOS_LAMBDA(const int radial_task, const int i_r) { const int i_theta = start_black_radials + radial_task * 2; - nodeApplyAscOrthoRadialTake(i_r, i_theta, grid, DirBC_Interior, x, rhs, temp, arr, att, art); + nodeApplyAscOrthoRadialTakeOuterBoundary(i_r, i_theta, grid, x, rhs, temp, arr, att, art); }); Kokkos::fence(); @@ -727,12 +492,12 @@ void ExtrapolatedSmootherTake::applyAscOrthoWhiteRadialSection(C ConstVector rhs, Vector temp) { - using extrapolated_smoother_take::nodeApplyAscOrthoRadialTake; + using extrapolated_smoother_take::nodeApplyAscOrthoRadialTakeInnerBoundary; using extrapolated_smoother_take::nodeApplyAscOrthoRadialTakeInterior; + using extrapolated_smoother_take::nodeApplyAscOrthoRadialTakeOuterBoundary; const PolarGrid& grid = ExtrapolatedSmoother::grid_; const LevelCacheType& level_cache = ExtrapolatedSmoother::level_cache_; - const bool DirBC_Interior = ExtrapolatedSmoother::DirBC_Interior_; assert(level_cache.cacheDomainGeometry()); @@ -753,7 +518,7 @@ void ExtrapolatedSmootherTake::applyAscOrthoWhiteRadialSection(C // Kokkos lambda function to execute for each point in the index space KOKKOS_LAMBDA(const int radial_task, const int i_r) { const int i_theta = start_white_radials + radial_task * 2; - nodeApplyAscOrthoRadialTake(i_r, i_theta, grid, DirBC_Interior, x, rhs, temp, arr, att, art); + nodeApplyAscOrthoRadialTakeInnerBoundary(i_r, i_theta, grid, x, rhs, temp, arr, att, art); }); Kokkos::parallel_for( @@ -765,7 +530,7 @@ void ExtrapolatedSmootherTake::applyAscOrthoWhiteRadialSection(C // Kokkos lambda function to execute for each point in the index space KOKKOS_LAMBDA(const int radial_task, const int i_r) { const int i_theta = start_white_radials + radial_task * 2; - nodeApplyAscOrthoRadialTakeInterior(i_r, i_theta, grid, DirBC_Interior, x, rhs, temp, arr, att, art); + nodeApplyAscOrthoRadialTakeInterior(i_r, i_theta, grid, x, rhs, temp, arr, att, art); }); Kokkos::parallel_for( @@ -777,7 +542,7 @@ void ExtrapolatedSmootherTake::applyAscOrthoWhiteRadialSection(C // Kokkos lambda function to execute for each point in the index space KOKKOS_LAMBDA(const int radial_task, const int i_r) { const int i_theta = start_white_radials + radial_task * 2; - nodeApplyAscOrthoRadialTake(i_r, i_theta, grid, DirBC_Interior, x, rhs, temp, arr, att, art); + nodeApplyAscOrthoRadialTakeOuterBoundary(i_r, i_theta, grid, x, rhs, temp, arr, att, art); }); Kokkos::fence(); diff --git a/include/Residual/ResidualGive/residualGive.h b/include/Residual/ResidualGive/residualGive.h index cb51ad8f..595c3b9b 100644 --- a/include/Residual/ResidualGive/residualGive.h +++ b/include/Residual/ResidualGive/residualGive.h @@ -10,7 +10,6 @@ class ResidualGive : public Residual { public: explicit ResidualGive(const PolarGrid& grid, const LevelCacheType& level_cache, const bool DirBC_Interior); - ~ResidualGive() override = default; void applySystemOperator(Vector result, ConstVector x) const final; void computeResidual(Vector result, ConstVector rhs, ConstVector x) const final; diff --git a/include/Residual/ResidualTake/applyATake.inl b/include/Residual/ResidualTake/applyATake.inl index d6f589fc..6eb47a29 100644 --- a/include/Residual/ResidualTake/applyATake.inl +++ b/include/Residual/ResidualTake/applyATake.inl @@ -3,12 +3,10 @@ namespace residual_take { -static KOKKOS_INLINE_FUNCTION void applySystemMatrixTakeInterior(const int i_r, const int i_theta, - const PolarGrid& grid, bool DirBC_Interior, - Vector& result, ConstVector& x, - ConstVector& arr, ConstVector& att, - ConstVector& art, ConstVector& detDF, - ConstVector& coeff_beta) +static KOKKOS_INLINE_FUNCTION void +applySystemMatrixTakeInterior(const int i_r, const int i_theta, const PolarGrid& grid, Vector& result, + ConstVector& x, ConstVector& arr, ConstVector& att, + ConstVector& art, ConstVector& detDF, ConstVector& coeff_beta) { KOKKOS_ASSERT(0 < i_r && i_r < grid.nr() - 1); @@ -138,6 +136,7 @@ void ResidualTake::applySystemOperator(Vector result, Co applySystemMatrixTakeBoundary(i_r, i_theta, grid, DirBC_Interior, result, x, arr, att, art, detDF, coeff_beta); }); + Kokkos::parallel_for( "Residual Take: Apply System Operator Interior (Circular)", Kokkos::MDRangePolicy>( // Rank of the index space @@ -146,8 +145,7 @@ void ResidualTake::applySystemOperator(Vector result, Co ), // Kokkos lambda function to execute for each point in the index space KOKKOS_LAMBDA(const int i_r, const int i_theta) { - applySystemMatrixTakeInterior(i_r, i_theta, grid, DirBC_Interior, result, x, arr, att, art, detDF, - coeff_beta); + applySystemMatrixTakeInterior(i_r, i_theta, grid, result, x, arr, att, art, detDF, coeff_beta); }); /* For loop matches radial access pattern */ @@ -159,9 +157,9 @@ void ResidualTake::applySystemOperator(Vector result, Co ), // Kokkos lambda function to execute for each point in the index space KOKKOS_LAMBDA(const int i_theta, const int i_r) { - applySystemMatrixTakeInterior(i_r, i_theta, grid, DirBC_Interior, result, x, arr, att, art, detDF, - coeff_beta); + applySystemMatrixTakeInterior(i_r, i_theta, grid, result, x, arr, att, art, detDF, coeff_beta); }); + Kokkos::parallel_for( "Residual Take: Apply System Operator Boundary (Radial)", Kokkos::MDRangePolicy>( // Rank of the index space diff --git a/include/Residual/ResidualTake/residualTake.h b/include/Residual/ResidualTake/residualTake.h index c6577590..baad5398 100644 --- a/include/Residual/ResidualTake/residualTake.h +++ b/include/Residual/ResidualTake/residualTake.h @@ -10,8 +10,6 @@ class ResidualTake : public Residual { public: explicit ResidualTake(const PolarGrid& grid, const LevelCacheType& level_cache, const bool DirBC_Interior); - KOKKOS_DEFAULTED_FUNCTION ResidualTake(const ResidualTake&) = default; - KOKKOS_DEFAULTED_FUNCTION ~ResidualTake() override = default; void applySystemOperator(Vector result, ConstVector x) const final; void computeResidual(Vector result, ConstVector rhs, ConstVector x) const final; diff --git a/include/Smoother/SmootherTake/applyAscOrtho.inl b/include/Smoother/SmootherTake/applyAscOrtho.inl index 41bcb986..df7fa8ea 100644 --- a/include/Smoother/SmootherTake/applyAscOrtho.inl +++ b/include/Smoother/SmootherTake/applyAscOrtho.inl @@ -4,11 +4,11 @@ namespace smoother_take { static KOKKOS_INLINE_FUNCTION void -nodeApplyAscOrthoCircleTakeInterior(const int i_r, const int i_theta, const PolarGrid& grid, bool DirBC_Interior, - ConstVector& x, ConstVector& rhs, Vector& result, - ConstVector& arr, ConstVector& att, ConstVector& art) +nodeApplyAscOrthoCircleTakeInterior(const int i_r, const int i_theta, const PolarGrid& grid, ConstVector& x, + ConstVector& rhs, Vector& result, ConstVector& arr, + ConstVector& att, ConstVector& art) { - KOKKOS_ASSERT(i_r >= 1 && i_r < grid.numberSmootherCircles()); + KOKKOS_ASSERT(0 < i_r && i_r < grid.numberSmootherCircles()); const int i_theta_M1 = grid.wrapThetaIndex(i_theta - 1); const int i_theta_P1 = grid.wrapThetaIndex(i_theta + 1); @@ -77,11 +77,11 @@ nodeApplyAscOrthoCircleTakeBoundary(const int i_r, const int i_theta, const Pola } static KOKKOS_INLINE_FUNCTION void -nodeApplyAscOrthoRadialTakeInterior(const int i_r, const int i_theta, const PolarGrid& grid, bool DirBC_Interior, - ConstVector& x, ConstVector& rhs, Vector& result, - ConstVector& arr, ConstVector& att, ConstVector& art) +nodeApplyAscOrthoRadialTakeInterior(const int i_r, const int i_theta, const PolarGrid& grid, ConstVector& x, + ConstVector& rhs, Vector& result, ConstVector& arr, + ConstVector& att, ConstVector& art) { - KOKKOS_ASSERT(i_r > grid.numberSmootherCircles() && i_r < grid.nr() - 2); + KOKKOS_ASSERT(grid.numberSmootherCircles() < i_r && 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); @@ -115,9 +115,9 @@ nodeApplyAscOrthoRadialTakeInterior(const int i_r, const int i_theta, const Pola } static KOKKOS_INLINE_FUNCTION void -nodeApplyAscOrthoRadialTakeBoundary(const int i_r, const int i_theta, const PolarGrid& grid, bool DirBC_Interior, - ConstVector& x, ConstVector& rhs, Vector& result, - ConstVector& arr, ConstVector& att, ConstVector& art) +nodeApplyAscOrthoRadialTakeBoundary(const int i_r, const int i_theta, const PolarGrid& grid, ConstVector& x, + ConstVector& rhs, Vector& result, ConstVector& arr, + ConstVector& att, ConstVector& art) { KOKKOS_ASSERT(i_r == grid.numberSmootherCircles() || i_r == grid.nr() - 2 || i_r == grid.nr() - 1); @@ -217,7 +217,7 @@ void SmootherTake::applyAscOrthoBlackCircleSection(ConstVector::applyAscOrthoBlackCircleSection(ConstVector::applyAscOrthoWhiteCircleSection(ConstVector::applyAscOrthoWhiteCircleSection(ConstVector::applyAscOrthoBlackRadialSection(ConstVector::applyAscOrthoBlackRadialSection(ConstVector::applyAscOrthoBlackRadialSection(ConstVector::applyAscOrthoWhiteRadialSection(ConstVector::applyAscOrthoWhiteRadialSection(ConstVector::applyAscOrthoWhiteRadialSection(ConstVector // Builds the A_sc smoother matrices and prepares the solvers. explicit SmootherTake(const PolarGrid& grid, const LevelCacheType& level_cache, bool DirBC_Interior); - KOKKOS_DEFAULTED_FUNCTION SmootherTake(const SmootherTake&) = default; - // Performs one full coupled smoothing sweep: // BC -> WC -> BR -> WR // using temp as RHS workspace. From 7192176d372f41dc4fb621c42d8713813fcb1c62 Mon Sep 17 00:00:00 2001 From: julianlitz Date: Sat, 22 Aug 2026 22:06:09 +0200 Subject: [PATCH 6/7] Fix --- scripts/tutorial/run.sh | 18 +++++++++--------- 1 file changed, 9 insertions(+), 9 deletions(-) diff --git a/scripts/tutorial/run.sh b/scripts/tutorial/run.sh index fb9ea8a5..1884a8f7 100755 --- a/scripts/tutorial/run.sh +++ b/scripts/tutorial/run.sh @@ -40,10 +40,10 @@ fi # Finest grid parameters R0=1e-8 Rmax=1.3 -nr_exp=4 -ntheta_exp=-1 -anisotropic_factor=3 -divideBy2=3 +nr_exp=5 +ntheta_exp=5 +anisotropic_factor=0 +divideBy2=2 # Interior boundary condition: # 0: Across-origin @@ -69,7 +69,7 @@ FMG_cycle=2 # V-Cycle(0), W-Cycle(1), F-Cycle(2) # Preconditioned Conjugate Gradient Method: # 0: GMGPolar as iterative solver # 1: GMGPolar solver as preconditioner for Conjugate Gradient (recommended) -PCG=1 +PCG=0 # Initial approximation for PCG: # 0: Initial approximation is set to residual -> no preconditioning # 1: FMG-approximation as initial guess (recommended) @@ -77,7 +77,7 @@ PCG_FMG=1 PCG_FMG_iterations=1 PCG_FMG_cycle=0 # V-Cycle(0), W-Cycle(1), F-Cycle(2) # Additional multigrid iterations after initial approximation to solve the linear system in each PCG iteration -PCG_MG_iterations=1 +PCG_MG_iterations=2 PCG_MG_cycle=0 # V-Cycle(0), W-Cycle(1), F-Cycle(2) # Extrapolation Method: @@ -98,10 +98,10 @@ postSmoothingSteps=1 multigridCycle=0 # Convergence criteria: -maxIterations=150 +maxIterations=1000 residualNormType=0 # L2-Norm(0) = 0, Weighted L2-Norm(1), Infinity-Norm(2) -absoluteTolerance=1e-10 -relativeTolerance=1e-10 +absoluteTolerance=-1e-12 +relativeTolerance=-1e-08 # Define additional geometry parameters kappa_eps=0.0 From e4b85738a17d5971f3188107991f9567e8101700 Mon Sep 17 00:00:00 2001 From: julianlitz Date: Sat, 22 Aug 2026 22:07:25 +0200 Subject: [PATCH 7/7] Fix --- .../applyAscOrtho.inl | 124 +++++++++++------- scripts/tutorial/run.sh | 20 +-- 2 files changed, 87 insertions(+), 57 deletions(-) diff --git a/include/ExtrapolatedSmoother/ExtrapolatedSmootherTake/applyAscOrtho.inl b/include/ExtrapolatedSmoother/ExtrapolatedSmootherTake/applyAscOrtho.inl index 1481bbfc..da709fc4 100644 --- a/include/ExtrapolatedSmoother/ExtrapolatedSmootherTake/applyAscOrtho.inl +++ b/include/ExtrapolatedSmoother/ExtrapolatedSmootherTake/applyAscOrtho.inl @@ -8,7 +8,7 @@ nodeApplyAscOrthoCircleTakeInterior(const int i_r, const int i_theta, const Pola ConstVector& rhs, Vector& result, ConstVector& arr, ConstVector& att, ConstVector& art) { - KOKKOS_ASSERT(0 < i_r && i_r <= grid.numberSmootherCircles()); + KOKKOS_ASSERT(0 < i_r && i_r < grid.numberSmootherCircles()); if (!(i_r & 1) && !(i_theta & 1)) { const int center = grid.index(i_r, i_theta); @@ -350,29 +350,44 @@ void ExtrapolatedSmootherTake::applyAscOrthoBlackCircleSection(C const int start_black_circles = (grid.numberSmootherCircles() % 2 == 0) ? 1 : 0; const int num_black_circles = (grid.numberSmootherCircles() - start_black_circles + 1) / 2; - Kokkos::parallel_for( - "ExtrapolatedSmootherTake: ApplyAscOrtho (Black Circular)", - Kokkos::MDRangePolicy>( // Rank of the index space - {0, 0}, // Starting point of the index space - {1, grid.ntheta()} // Ending point of the index space - ), - // Kokkos lambda function to execute for each point in the index space - KOKKOS_LAMBDA(const int circle_task, const int i_theta) { - int i_r = start_black_circles + circle_task * 2; - nodeApplyAscOrthoCircleTakeBoundary(i_r, i_theta, grid, DirBC_Interior, x, rhs, temp, arr, att, art); - }); - - Kokkos::parallel_for( - "ExtrapolatedSmootherTake: ApplyAscOrtho (Black Circular)", - Kokkos::MDRangePolicy>( // Rank of the index space - {1, 0}, // Starting point of the index space - {num_black_circles, grid.ntheta()} // Ending point of the index space - ), - // Kokkos lambda function to execute for each point in the index space - KOKKOS_LAMBDA(const int circle_task, const int i_theta) { - int i_r = start_black_circles + circle_task * 2; - nodeApplyAscOrthoCircleTakeInterior(i_r, i_theta, grid, x, rhs, temp, arr, att, art); - }); + if (start_black_circles == 0) { + Kokkos::parallel_for( + "ExtrapolatedSmootherTake: ApplyAscOrtho (Black Circular)", + Kokkos::MDRangePolicy>( // Rank of the index space + {0, 0}, // Starting point of the index space + {1, grid.ntheta()} // Ending point of the index space + ), + // Kokkos lambda function to execute for each point in the index space + KOKKOS_LAMBDA(const int circle_task, const int i_theta) { + int i_r = start_black_circles + circle_task * 2; + nodeApplyAscOrthoCircleTakeBoundary(i_r, i_theta, grid, DirBC_Interior, x, rhs, temp, arr, att, art); + }); + + Kokkos::parallel_for( + "ExtrapolatedSmootherTake: ApplyAscOrtho (Black Circular)", + Kokkos::MDRangePolicy>( // Rank of the index space + {1, 0}, // Starting point of the index space + {num_black_circles, grid.ntheta()} // Ending point of the index space + ), + // Kokkos lambda function to execute for each point in the index space + KOKKOS_LAMBDA(const int circle_task, const int i_theta) { + int i_r = start_black_circles + circle_task * 2; + nodeApplyAscOrthoCircleTakeInterior(i_r, i_theta, grid, x, rhs, temp, arr, att, art); + }); + } + else { + Kokkos::parallel_for( + "ExtrapolatedSmootherTake: ApplyAscOrtho (Black Circular)", + Kokkos::MDRangePolicy>( // Rank of the index space + {0, 0}, // Starting point of the index space + {num_black_circles, grid.ntheta()} // Ending point of the index space + ), + // Kokkos lambda function to execute for each point in the index space + KOKKOS_LAMBDA(const int circle_task, const int i_theta) { + int i_r = start_black_circles + circle_task * 2; + nodeApplyAscOrthoCircleTakeInterior(i_r, i_theta, grid, x, rhs, temp, arr, att, art); + }); + } Kokkos::fence(); } @@ -399,29 +414,44 @@ void ExtrapolatedSmootherTake::applyAscOrthoWhiteCircleSection(C const int start_white_circles = (grid.numberSmootherCircles() % 2 == 0) ? 0 : 1; const int num_white_circles = (grid.numberSmootherCircles() - start_white_circles + 1) / 2; - Kokkos::parallel_for( - "ExtrapolatedSmootherTake: ApplyAscOrtho (White Circular)", - Kokkos::MDRangePolicy>( // Rank of the index space - {0, 0}, // Starting point of the index space - {1, grid.ntheta()} // Ending point of the index space - ), - // Kokkos lambda function to execute for each point in the index space - KOKKOS_LAMBDA(const int circle_task, const int i_theta) { - const int i_r = start_white_circles + circle_task * 2; - nodeApplyAscOrthoCircleTakeBoundary(i_r, i_theta, grid, DirBC_Interior, x, rhs, temp, arr, att, art); - }); - - Kokkos::parallel_for( - "ExtrapolatedSmootherTake: ApplyAscOrtho (White Circular)", - Kokkos::MDRangePolicy>( // Rank of the index space - {1, 0}, // Starting point of the index space - {num_white_circles, grid.ntheta()} // Ending point of the index space - ), - // Kokkos lambda function to execute for each point in the index space - KOKKOS_LAMBDA(const int circle_task, const int i_theta) { - const int i_r = start_white_circles + circle_task * 2; - nodeApplyAscOrthoCircleTakeInterior(i_r, i_theta, grid, x, rhs, temp, arr, att, art); - }); + if (start_white_circles == 0) { + Kokkos::parallel_for( + "ExtrapolatedSmootherTake: ApplyAscOrtho (White Circular)", + Kokkos::MDRangePolicy>( // Rank of the index space + {0, 0}, // Starting point of the index space + {1, grid.ntheta()} // Ending point of the index space + ), + // Kokkos lambda function to execute for each point in the index space + KOKKOS_LAMBDA(const int circle_task, const int i_theta) { + const int i_r = start_white_circles + circle_task * 2; + nodeApplyAscOrthoCircleTakeBoundary(i_r, i_theta, grid, DirBC_Interior, x, rhs, temp, arr, att, art); + }); + + Kokkos::parallel_for( + "ExtrapolatedSmootherTake: ApplyAscOrtho (White Circular)", + Kokkos::MDRangePolicy>( // Rank of the index space + {1, 0}, // Starting point of the index space + {num_white_circles, grid.ntheta()} // Ending point of the index space + ), + // Kokkos lambda function to execute for each point in the index space + KOKKOS_LAMBDA(const int circle_task, const int i_theta) { + const int i_r = start_white_circles + circle_task * 2; + nodeApplyAscOrthoCircleTakeInterior(i_r, i_theta, grid, x, rhs, temp, arr, att, art); + }); + } + else { + Kokkos::parallel_for( + "ExtrapolatedSmootherTake: ApplyAscOrtho (White Circular)", + Kokkos::MDRangePolicy>( // Rank of the index space + {0, 0}, // Starting point of the index space + {num_white_circles, grid.ntheta()} // Ending point of the index space + ), + // Kokkos lambda function to execute for each point in the index space + KOKKOS_LAMBDA(const int circle_task, const int i_theta) { + const int i_r = start_white_circles + circle_task * 2; + nodeApplyAscOrthoCircleTakeInterior(i_r, i_theta, grid, x, rhs, temp, arr, att, art); + }); + } Kokkos::fence(); } diff --git a/scripts/tutorial/run.sh b/scripts/tutorial/run.sh index 1884a8f7..572d4dd4 100755 --- a/scripts/tutorial/run.sh +++ b/scripts/tutorial/run.sh @@ -40,10 +40,10 @@ fi # Finest grid parameters R0=1e-8 Rmax=1.3 -nr_exp=5 -ntheta_exp=5 -anisotropic_factor=0 -divideBy2=2 +nr_exp=4 +ntheta_exp=-1 +anisotropic_factor=3 +divideBy2=3 # Interior boundary condition: # 0: Across-origin @@ -69,7 +69,7 @@ FMG_cycle=2 # V-Cycle(0), W-Cycle(1), F-Cycle(2) # Preconditioned Conjugate Gradient Method: # 0: GMGPolar as iterative solver # 1: GMGPolar solver as preconditioner for Conjugate Gradient (recommended) -PCG=0 +PCG=1 # Initial approximation for PCG: # 0: Initial approximation is set to residual -> no preconditioning # 1: FMG-approximation as initial guess (recommended) @@ -77,7 +77,7 @@ PCG_FMG=1 PCG_FMG_iterations=1 PCG_FMG_cycle=0 # V-Cycle(0), W-Cycle(1), F-Cycle(2) # Additional multigrid iterations after initial approximation to solve the linear system in each PCG iteration -PCG_MG_iterations=2 +PCG_MG_iterations=1 PCG_MG_cycle=0 # V-Cycle(0), W-Cycle(1), F-Cycle(2) # Extrapolation Method: @@ -98,10 +98,10 @@ postSmoothingSteps=1 multigridCycle=0 # Convergence criteria: -maxIterations=1000 +maxIterations=150 residualNormType=0 # L2-Norm(0) = 0, Weighted L2-Norm(1), Infinity-Norm(2) -absoluteTolerance=-1e-12 -relativeTolerance=-1e-08 +absoluteTolerance=1e-10 +relativeTolerance=1e-10 # Define additional geometry parameters kappa_eps=0.0 @@ -168,4 +168,4 @@ export OMP_NUM_THREADS=$maxOpenMPThreads --maxIterations $maxIterations \ --residualNormType $residualNormType \ --absoluteTolerance $absoluteTolerance \ - --relativeTolerance $relativeTolerance + --relativeTolerance $relativeTolerance \ No newline at end of file