From ef2a35c394087bf63ca370c001d46774114ac983 Mon Sep 17 00:00:00 2001 From: julianlitz Date: Mon, 17 Aug 2026 00:14:46 +0200 Subject: [PATCH 1/4] 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/4] 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/4] 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/4] 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);