From 8c132f1bd30f973ebe7a8f026f7012bcc6f37aa1 Mon Sep 17 00:00:00 2001 From: julianlitz Date: Wed, 13 May 2026 19:39:51 +0200 Subject: [PATCH 01/10] Test --- .../Smoother/SmootherGive/applyAscOrtho.inl | 850 ++++++++++-------- include/Smoother/SmootherGive/smootherGive.h | 10 +- .../Smoother/SmootherGive/smootherGive.inl | 105 +-- 3 files changed, 483 insertions(+), 482 deletions(-) diff --git a/include/Smoother/SmootherGive/applyAscOrtho.inl b/include/Smoother/SmootherGive/applyAscOrtho.inl index f91f5366e..a8a890394 100644 --- a/include/Smoother/SmootherGive/applyAscOrtho.inl +++ b/include/Smoother/SmootherGive/applyAscOrtho.inl @@ -3,71 +3,56 @@ namespace smoother_give { -static inline void nodeApplyAscOrthoCircleGive(int i_r, int i_theta, const PolarGrid& grid, bool DirBC_Interior, - SmootherColor smoother_color, ConstVector& x, - ConstVector& rhs, Vector& result, double arr, double att, - double art, double detDF, double coeff_beta) +template +static inline void nodeApplyAscOrthoCircleGiveInside(int i_r, int i_theta, const PolarGrid& grid, + const LevelCacheType& level_cache, bool DirBC_Interior, + ConstVector& x, ConstVector& rhs, + Vector& result) { - assert(i_r >= 0 && i_r <= grid.numberSmootherCircles()); + assert(i_r >= 0 && i_r < grid.numberSmootherCircles()); - bool isOddNumberSmootherCircles = (grid.numberSmootherCircles() & 1); - bool isOddRadialIndex = (i_r & 1); + /* ---------------------------------------- */ + /* Compute or retrieve stencil coefficients */ + /* ---------------------------------------- */ + const int center = grid.index(i_r, i_theta); + const double radius = grid.radius(i_r); + const double theta = grid.theta(i_theta); - SmootherColor node_color = - (isOddNumberSmootherCircles == isOddRadialIndex) ? SmootherColor::White : SmootherColor::Black; + double coeff_beta, arr, att, art, detDF; + level_cache.obtainValues(i_r, i_theta, center, radius, theta, coeff_beta, arr, att, art, detDF); /* -------------------- */ /* Node in the interior */ /* -------------------- */ if (i_r > 0 && i_r < grid.numberSmootherCircles()) { - double h1 = grid.radialSpacing(i_r - 1); - double h2 = grid.radialSpacing(i_r); - double k1 = grid.angularSpacing(i_theta - 1); - double k2 = grid.angularSpacing(i_theta); - - double coeff1 = 0.5 * (k1 + k2) / h1; - double coeff2 = 0.5 * (k1 + k2) / h2; - double coeff3 = 0.5 * (h1 + h2) / k1; - double coeff4 = 0.5 * (h1 + h2) / k2; - - int i_theta_M1 = grid.wrapThetaIndex(i_theta - 1); - int i_theta_P1 = grid.wrapThetaIndex(i_theta + 1); - - int center = grid.index(i_r, i_theta); - int left = grid.index(i_r - 1, i_theta); - int right = grid.index(i_r + 1, i_theta); - int bottom = grid.index(i_r, i_theta_M1); - int top = grid.index(i_r, i_theta_P1); - - /* -------------------- */ - /* Inside Section Parts */ - if (node_color == smoother_color) { - /* Fill result(i,j) */ - result[center] -= (-coeff1 * arr * x[left] /* Left */ - - coeff2 * arr * x[right]); /* Right */ - /* Fill result(i,j-1) */ - result[bottom] -= (-0.25 * art * x[right] /* Top Right */ - + 0.25 * art * x[left]); /* Top Left */ - /* Fill result(i,j+1) */ - result[top] -= (+0.25 * art * x[right] /* Bottom Right */ - - 0.25 * art * x[left]); /* Bottom Left */ - } - /* --------------------- */ - /* Outside Section Parts */ - else if (node_color != smoother_color) { - /* Fill result(i-1,j) */ - if (i_r > 1 || !DirBC_Interior) { - result[left] -= (-coeff1 * arr * x[center] /* Right */ - - 0.25 * art * x[top] /* Top Right */ - + 0.25 * art * x[bottom]); /* Bottom Right */ - } - /* Fill result(i+1,j) */ - if (i_r < grid.numberSmootherCircles() - 1) { - result[right] -= (-coeff2 * arr * x[center] /* Left */ - + 0.25 * art * x[top] /* Top Left */ - - 0.25 * art * x[bottom]); /* Bottom Left */ - } - } + 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 center = grid.index(i_r, i_theta); + const int left = grid.index(i_r - 1, i_theta); + const int right = grid.index(i_r + 1, i_theta); + const int bottom = grid.index(i_r, i_theta_M1); + const int top = grid.index(i_r, i_theta_P1); + + /* Fill result(i,j) */ + result[center] -= (-coeff1 * arr * x[left] /* Left */ + - coeff2 * arr * x[right]); /* Right */ + /* Fill result(i,j-1) */ + result[bottom] -= (-0.25 * art * x[right] /* Top Right */ + + 0.25 * art * x[left]); /* Top Left */ + /* Fill result(i,j+1) */ + result[top] -= (+0.25 * art * x[right] /* Bottom Right */ + - 0.25 * art * x[left]); /* Bottom Left */ } /* -------------------- */ /* Node on the boundary */ @@ -77,215 +62,183 @@ static inline void nodeApplyAscOrthoCircleGive(int i_r, int i_theta, const Polar /* Case 1: Dirichlet boundary on the inner boundary */ /* ------------------------------------------------ */ if (DirBC_Interior) { - double h2 = grid.radialSpacing(i_r); - double k1 = grid.angularSpacing(i_theta - 1); - double k2 = grid.angularSpacing(i_theta); - - double coeff2 = 0.5 * (k1 + k2) / h2; - - int i_theta_M1 = grid.wrapThetaIndex(i_theta - 1); - int i_theta_P1 = grid.wrapThetaIndex(i_theta + 1); - - int center = grid.index(i_r, i_theta); - int right = grid.index(i_r + 1, i_theta); - int bottom = grid.index(i_r, i_theta_M1); - int top = grid.index(i_r, i_theta_P1); - - /* -------------------- */ - /* Inside Section Parts */ - if (node_color == smoother_color) { - /* Nothing to be done here */ - } - /* --------------------- */ - /* Outside Section Parts */ - else if (node_color != smoother_color) { - /* Fill result(i+1,j) */ - result[right] -= (-coeff2 * arr * x[center] /* Left */ - + 0.25 * art * x[top] /* Top Left */ - - 0.25 * art * x[bottom]); /* Bottom Left */ - } + /* Nothing to be done here */ } 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. - double h1 = 2.0 * grid.radius(0); - double h2 = grid.radialSpacing(i_r); - double k1 = grid.angularSpacing(i_theta - 1); - double k2 = grid.angularSpacing(i_theta); - - double coeff1 = 0.5 * (k1 + k2) / h1; - double coeff2 = 0.5 * (k1 + k2) / h2; - double coeff3 = 0.5 * (h1 + h2) / k1; - 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 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 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 right = grid.index(i_r + 1, i_theta); - /* -------------------- */ - /* Inside Section Parts */ - if (node_color == smoother_color) { - /* Fill result(i,j) */ - result[center] -= ( - /* - coeff1 * arr * x[left] // Left: Not in Asc_ortho */ - -coeff2 * arr * x[right]); /* Right */ - /* Fill result(i,j-1) */ - result[bottom] -= (-0.25 * art * x[right]); /* Top Right */ - /* + 0.25 * art * x[left]; // Top Left: REMOVED DUE TO ARTIFICAL 7 POINT STENCIL */ - /* Fill result(i,j+1) */ - result[top] -= (+0.25 * art * x[right]); /* Bottom Right */ - /* - 0.25 * art * x[left]; // Bottom Left: REMOVED DUE TO ARTIFICAL 7 POINT STENCIL */ - } - /* --------------------- */ - /* Outside Section Parts */ - else if (node_color != smoother_color) { - /* Fill result(i+1,j) */ - result[right] -= (-coeff2 * arr * x[center] /* Left */ - + 0.25 * art * x[top] /* Top Left */ - - 0.25 * art * x[bottom]); /* Bottom Left */ - } + /* Fill result(i,j) */ + result[center] -= (-coeff2 * arr * x[right]); /* Right */ + + /* Fill result(i,j-1) */ + result[bottom] -= (-0.25 * art * x[right]); /* Top Right */ + + /* Fill result(i,j+1) */ + result[top] -= (+0.25 * art * x[right]); /* Bottom Right */ } } - /* ----------------------------- */ - /* Node next to circular section */ - /* ----------------------------- */ - else if (i_r == grid.numberSmootherCircles()) { - assert(node_color == SmootherColor::White); - - if (smoother_color == SmootherColor::Black) { - double h1 = grid.radialSpacing(i_r - 1); - double k1 = grid.angularSpacing(i_theta - 1); - double k2 = grid.angularSpacing(i_theta); - - double coeff1 = 0.5 * (k1 + k2) / h1; - - int i_theta_M1 = grid.wrapThetaIndex(i_theta - 1); - int i_theta_P1 = grid.wrapThetaIndex(i_theta + 1); - - int center = grid.index(i_r, i_theta); - int left = grid.index(i_r - 1, i_theta); - int bottom = grid.index(i_r, i_theta_M1); - int top = grid.index(i_r, i_theta_P1); - - /* --------------------- */ - /* Outside Section Parts */ - /* Fill result(i-1,j) */ +} + +template +static inline void nodeApplyAscOrthoCircleGiveOutside(int i_r, int i_theta, const PolarGrid& grid, + const LevelCacheType& level_cache, bool DirBC_Interior, + ConstVector& x, ConstVector& rhs, + Vector& result) +{ + assert(0 <= i_r && i_r <= grid.numberSmootherCircles()); + + /* ---------------------------------------- */ + /* Compute or retrieve stencil coefficients */ + /* ---------------------------------------- */ + const int center = grid.index(i_r, i_theta); + const double radius = grid.radius(i_r); + const double theta = grid.theta(i_theta); + + double coeff_beta, arr, att, art, detDF; + level_cache.obtainValues(i_r, i_theta, center, radius, theta, coeff_beta, arr, att, art, detDF); + + if (0 <= i_r && 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 left = grid.index(i_r - 1, i_theta); + const int right = grid.index(i_r + 1, i_theta); + const int bottom = grid.index(i_r, i_theta_M1); + const int top = grid.index(i_r, i_theta_P1); + + /* Fill result(i-1,j) */ + if (i_r > 1 || !DirBC_Interior) { result[left] -= (-coeff1 * arr * x[center] /* Right */ - 0.25 * art * x[top] /* Top Right */ + 0.25 * art * x[bottom]); /* Bottom Right */ } + /* Fill result(i+1,j) */ + if (i_r < grid.numberSmootherCircles() - 1) { + result[right] -= (-coeff2 * arr * x[center] /* Left */ + + 0.25 * art * x[top] /* Top Left */ + - 0.25 * art * x[bottom]); /* Bottom Left */ + } } } -static inline void nodeApplyAscOrthoRadialGive(int i_r, int i_theta, const PolarGrid& grid, bool DirBC_Interior, - SmootherColor smoother_color, ConstVector& x, - ConstVector& rhs, Vector& result, double arr, double att, - double art, double detDF, double coeff_beta) +template +static inline void nodeApplyAscOrthoRadialGiveInside(int i_r, int i_theta, const PolarGrid& grid, + const LevelCacheType& level_cache, bool DirBC_Interior, + ConstVector& x, ConstVector& rhs, + Vector& result) { - assert(i_r >= grid.numberSmootherCircles() - 1 && i_r < grid.nr()); - - SmootherColor node_color = (i_theta & 1) ? SmootherColor::White : SmootherColor::Black; - - /* -------------------- */ - /* Node in the interior */ - /* -------------------- */ - if (i_r > grid.numberSmootherCircles() && i_r < grid.nr() - 2) { - double h1 = grid.radialSpacing(i_r - 1); - double h2 = grid.radialSpacing(i_r); - double k1 = grid.angularSpacing(i_theta - 1); - double k2 = grid.angularSpacing(i_theta); + assert(grid.numberSmootherCircles() - 1 <= i_r && i_r < grid.nr()); + + /* ---------------------------------------- */ + /* Compute or retrieve stencil coefficients */ + /* ---------------------------------------- */ + const int center = grid.index(i_r, i_theta); + const double radius = grid.radius(i_r); + const double theta = grid.theta(i_theta); + + double coeff_beta, arr, att, art, detDF; + level_cache.obtainValues(i_r, i_theta, center, radius, theta, coeff_beta, arr, att, art, detDF); + + /* Angular neighbours (wrapping always valid) */ + 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 top = grid.index(i_r, i_theta_P1); + const int left = grid.index(i_r - 1, i_theta); + + /* h1, k1, k2 are always valid at this point */ + const double h1 = grid.radialSpacing(i_r - 1); + const double k1 = grid.angularSpacing(i_theta_M1); + const double k2 = grid.angularSpacing(i_theta); + const double coeff1 = 0.5 * (k1 + k2) / h1; + + /* ------------------- */ + /* Fill result(center) */ + /* ------------------- */ + if (grid.numberSmootherCircles() <= i_r && i_r < grid.nr() - 1) { + const double h2 = grid.radialSpacing(i_r); // safe: i_r < nr-1 + const double coeff3 = 0.5 * (h1 + h2) / k1; + const double coeff4 = 0.5 * (h1 + h2) / k2; + result[center] -= (-coeff3 * att * x[bottom] /* Bottom */ + - coeff4 * att * x[top]); /* Top */ + } - double coeff1 = 0.5 * (k1 + k2) / h1; - double coeff2 = 0.5 * (k1 + k2) / h2; - double coeff3 = 0.5 * (h1 + h2) / k1; - double coeff4 = 0.5 * (h1 + h2) / k2; + /* ----------------- */ + /* Fill result(left) */ + /* ----------------- */ + if (grid.numberSmootherCircles() < i_r && i_r < grid.nr()) { + result[left] -= (-0.25 * art * x[top] /* Top Right */ + + 0.25 * art * x[bottom]); /* Bottom Right */ + } - int i_theta_M1 = grid.wrapThetaIndex(i_theta - 1); - int i_theta_P1 = grid.wrapThetaIndex(i_theta + 1); + /* ------------------- */ + /* Fill result(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 */ + } - int center = grid.index(i_r, i_theta); - int left = grid.index(i_r - 1, i_theta); - int right = grid.index(i_r + 1, i_theta); - int bottom = grid.index(i_r, i_theta_M1); - int top = grid.index(i_r, i_theta_P1); + /* --------------- */ + /* Symmetry shifts */ + /* --------------- */ + if (i_r == grid.nr() - 2) { + 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 * rhs[right]; /* Right: Symmetry shift! */ + } - /* -------------------- */ - /* Inside Section Parts */ - if (node_color == smoother_color) { - /* Fill result(i,j) */ - result[center] -= (-coeff3 * att * x[bottom] /* Bottom */ - - coeff4 * att * x[top] /* Top */ - ); /* Fill result(i-1,j) */ - result[left] -= (-0.25 * art * x[top] /* Top Right */ - + 0.25 * art * x[bottom]); /* Bottom Right */ - /* Fill result(i+1,j) */ - result[right] -= (+0.25 * art * x[top] /* Top Left */ - - 0.25 * art * x[bottom]); /* Bottom Left */ - } - /* --------------------- */ - /* Outside Section Parts */ - else if (node_color != smoother_color) { - /* 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 */ - /* 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 */ - } + if (i_r == grid.nr() - 1) { + result[left] -= -coeff1 * arr * rhs[center]; /* Right: Symmetry shift! */ } +} - /* -------------------------------- */ - /* Node left next to radial section */ - /* -------------------------------- */ - else if (i_r == grid.numberSmootherCircles() - 1) { - double h1 = grid.radialSpacing(i_r - 1); - double h2 = grid.radialSpacing(i_r); - double k1 = grid.angularSpacing(i_theta - 1); - double k2 = grid.angularSpacing(i_theta); +template +static inline void nodeApplyAscOrthoRadialGiveOutside(int i_r, int i_theta, const PolarGrid& grid, + const LevelCacheType& level_cache, bool DirBC_Interior, + ConstVector& x, ConstVector& rhs, + Vector& result) +{ + assert(grid.numberSmootherCircles() <= i_r && i_r < grid.nr()); - double coeff1 = 0.5 * (k1 + k2) / h1; - double coeff2 = 0.5 * (k1 + k2) / h2; - double coeff3 = 0.5 * (h1 + h2) / k1; - double coeff4 = 0.5 * (h1 + h2) / k2; + /* ---------------------------------------- */ + /* Compute or retrieve stencil coefficients */ + /* ---------------------------------------- */ + const int center = grid.index(i_r, i_theta); + const double radius = grid.radius(i_r); + const double theta = grid.theta(i_theta); - int i_theta_M1 = grid.wrapThetaIndex(i_theta - 1); - int i_theta_P1 = grid.wrapThetaIndex(i_theta + 1); + double coeff_beta, arr, att, art, detDF; + level_cache.obtainValues(i_r, i_theta, center, radius, theta, coeff_beta, arr, att, art, detDF); - int center = grid.index(i_r, i_theta); - int right = grid.index(i_r + 1, i_theta); - int bottom = grid.index(i_r, i_theta_M1); - int top = grid.index(i_r, i_theta_P1); - - /* -------------------- */ - /* Inside Section Parts */ - if (node_color == smoother_color) { - /* Fill result(i+1,j) */ - result[right] -= (-coeff2 * arr * x[center] /* Left */ - + 0.25 * art * x[top] /* Top Left */ - - 0.25 * art * x[bottom]); /* Bottom Left */ - } - /* --------------------- */ - /* Outside Section Parts */ - else if (node_color != smoother_color) { - /* Nothing to be done here */ - } - } - /* --------------------------- */ - /* Node next to circle section */ - /* --------------------------- */ - else if (i_r == grid.numberSmootherCircles()) { + if (grid.numberSmootherCircles() <= i_r && i_r < grid.nr() - 1) { double h1 = grid.radialSpacing(i_r - 1); double h2 = grid.radialSpacing(i_r); double k1 = grid.angularSpacing(i_theta - 1); @@ -299,184 +252,327 @@ static inline void nodeApplyAscOrthoRadialGive(int i_r, int i_theta, const Polar int i_theta_M1 = grid.wrapThetaIndex(i_theta - 1); int i_theta_P1 = grid.wrapThetaIndex(i_theta + 1); - int center = grid.index(i_r, i_theta); int left = grid.index(i_r - 1, i_theta); int right = grid.index(i_r + 1, i_theta); int bottom = grid.index(i_r, i_theta_M1); int top = grid.index(i_r, i_theta_P1); - /* -------------------- */ - /* Inside Section Parts */ - if (node_color == smoother_color) { - /* Fill result(i,j) */ - result[center] -= (-coeff1 * arr * x[left] /* Left */ - - coeff3 * att * x[bottom] /* Bottom */ - - coeff4 * att * x[top] /* Top */ - ); - /* Fill result(i+1,j) */ - result[right] -= (+0.25 * art * x[top] /* Top Left */ - - 0.25 * art * x[bottom]); /* Bottom Left */ - } - /* --------------------- */ - /* Outside Section Parts */ - else if (node_color != smoother_color) { - /* 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 */ - /* 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 */ - } + /* 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 */ + /* 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 */ } - /* ------------------------------- */ - /* Node next to dirichlet boundary */ - /* ------------------------------- */ - else if (i_r == grid.nr() - 2) { - double h1 = grid.radialSpacing(i_r - 1); - double h2 = grid.radialSpacing(i_r); - double k1 = grid.angularSpacing(i_theta - 1); - double k2 = grid.angularSpacing(i_theta); +} - double coeff1 = 0.5 * (k1 + k2) / h1; - double coeff2 = 0.5 * (k1 + k2) / h2; - double coeff3 = 0.5 * (h1 + h2) / k1; - double coeff4 = 0.5 * (h1 + h2) / k2; +} // namespace smoother_give - int i_theta_M1 = grid.wrapThetaIndex(i_theta - 1); - int i_theta_P1 = grid.wrapThetaIndex(i_theta + 1); +template +void SmootherGive::applyAscOrthoBlackCircleSection(ConstVector x, ConstVector rhs, + Vector temp) +{ + using smoother_give::nodeApplyAscOrthoCircleGiveInside; + using smoother_give::nodeApplyAscOrthoCircleGiveOutside; - int center = grid.index(i_r, i_theta); - int left = grid.index(i_r - 1, i_theta); - int right = grid.index(i_r + 1, i_theta); - int bottom = grid.index(i_r, i_theta_M1); - int top = grid.index(i_r, i_theta_P1); + auto getBatchCount = [](int start, int end, int offset) { + if (start >= end) { + return 0; + } + return (end - start + offset - 1) / offset; + }; - /* -------------------- */ - /* Inside Section Parts */ - if (node_color == smoother_color) { - /* Fill result(i,j) */ - result[center] -= (-coeff3 * att * x[bottom] /* Bottom */ - - coeff4 * att * x[top] /* Top */ - ); - /* Fill result(i-1,j) */ - result[left] -= (-0.25 * art * x[top] /* Top Right */ - + 0.25 * art * x[bottom]); /* Bottom Right */ + const PolarGrid& grid = Smoother::grid_; + const LevelCacheType& level_cache = Smoother::level_cache_; + const bool DirBC_Interior = Smoother::DirBC_Interior_; - // "Right" is part of the radial Asc smoother matrices, - // but is shifted over to the rhs to make the radial Asc smoother matrices symmetric. - // Note that the circle Asc smoother matrices are symmetric by default. - // Note that rhs[right] contains the correct boundary value of u_D. - result[center] -= -coeff2 * arr * rhs[right]; /* Right: Symmetry shift! */ - } - /* --------------------- */ - /* Outside Section Parts */ - else if (node_color != smoother_color) { - /* 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 */ - /* 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 */ - } + /* ----------------------------------------------- */ + /* 1. Black-Circle update (u_bc): */ + /* A_bc * u_bc = f_bc − A_bc^ortho * u_bc^ortho */ + /* ----------------------------------------------- */ + { + /* Inside Black Section */ + const int start = (grid.numberSmootherCircles() % 2 == 0) ? 1 : 0; + const int end = grid.numberSmootherCircles(); + const int offset = 2; + + Kokkos::parallel_for( + "SmootherGive: ApplyAscOrtho (Black Circle - Inside)", + Kokkos::RangePolicy(0, getBatchCount(start, end, offset)), + KOKKOS_LAMBDA(const int circle_task) { + const int i_r = start + circle_task * offset; + for (int i_theta = 0; i_theta < grid.ntheta(); i_theta++) { + nodeApplyAscOrthoCircleGiveInside(i_r, i_theta, grid, level_cache, DirBC_Interior, x, rhs, temp); + } + }); + Kokkos::fence(); } - /* ------------------------------ */ - /* Node on the dirichlet boundary */ - /* ------------------------------ */ - else if (i_r == grid.nr() - 1) { - double h1 = grid.radialSpacing(i_r - 1); - double k1 = grid.angularSpacing(i_theta - 1); - double k2 = grid.angularSpacing(i_theta); - double coeff1 = 0.5 * (k1 + k2) / h1; + { + /* Outside Black Section (Part 1)*/ + const int start = (grid.numberSmootherCircles() % 2 == 0) ? 0 : 1; + const int end = grid.numberSmootherCircles() + 1; + const int offset = 4; + Kokkos::parallel_for( + "SmootherGive: ApplyAscOrtho (Black Circle - Outside: Part 1)", + Kokkos::RangePolicy(0, getBatchCount(start, end, offset)), + KOKKOS_LAMBDA(const int circle_task) { + const int i_r = start + circle_task * offset; + for (int i_theta = 0; i_theta < grid.ntheta(); i_theta++) { + nodeApplyAscOrthoCircleGiveOutside(i_r, i_theta, grid, level_cache, DirBC_Interior, x, rhs, temp); + } + }); + Kokkos::fence(); + } - int i_theta_M1 = grid.wrapThetaIndex(i_theta - 1); - int i_theta_P1 = grid.wrapThetaIndex(i_theta + 1); + { + /* Outside Black Section (Part 2)*/ + const int start = (grid.numberSmootherCircles() % 2 == 0) ? 2 : 3; + const int end = grid.numberSmootherCircles() + 1; + const int offset = 4; + Kokkos::parallel_for( + "SmootherGive: ApplyAscOrtho (Black Circle - Outside: Part 2)", + Kokkos::RangePolicy(0, getBatchCount(start, end, offset)), + KOKKOS_LAMBDA(const int circle_task) { + const int i_r = start + circle_task * offset; + for (int i_theta = 0; i_theta < grid.ntheta(); i_theta++) { + nodeApplyAscOrthoCircleGiveOutside(i_r, i_theta, grid, level_cache, DirBC_Interior, x, rhs, temp); + } + }); + Kokkos::fence(); + } +} - int center = grid.index(i_r, i_theta); - int left = grid.index(i_r - 1, i_theta); - int bottom = grid.index(i_r, i_theta_M1); - int top = grid.index(i_r, i_theta_P1); +template +void SmootherGive::applyAscOrthoWhiteCircleSection(ConstVector x, ConstVector rhs, + Vector temp) +{ + using smoother_give::nodeApplyAscOrthoCircleGiveInside; + using smoother_give::nodeApplyAscOrthoCircleGiveOutside; - /* -------------------- */ - /* Inside Section Parts */ - if (node_color == smoother_color) { - /* Fill result(i-1,j) */ - result[left] -= (-0.25 * art * x[top] /* Top Right */ - + 0.25 * art * x[bottom] /* Bottom Right */ - ); - // "Right" is part of the radial Asc smoother matrices, - // but is shifted over to the rhs to make the radial Asc smoother matrices symmetric. - // Note that the circle Asc smoother matrices are symmetric by default. - // Note that rhs[center] contains the correct boundary value of u_D. - result[left] -= (-coeff1 * arr * rhs[center] /* Right */ - ); - } - /* --------------------- */ - /* Outside Section Parts */ - else if (node_color != smoother_color) { - /* Nothing to be done here */ + auto getBatchCount = [](int start, int end, int offset) { + if (start >= end) { + return 0; } + return (end - start + offset - 1) / offset; + }; + + const PolarGrid& grid = Smoother::grid_; + const LevelCacheType& level_cache = Smoother::level_cache_; + const bool DirBC_Interior = Smoother::DirBC_Interior_; + + /* ----------------------------------------------- */ + /* 2. White-Circle update (u_wc): */ + /* A_wc * u_wc = f_wc − A_wc^ortho * u_wc^ortho */ + /* ----------------------------------------------- */ + + { + /* Inside White Section */ + const int start = (grid.numberSmootherCircles() % 2 == 0) ? 0 : 1; + const int end = grid.numberSmootherCircles(); + const int offset = 2; + + Kokkos::parallel_for( + "SmootherGive: ApplyAscOrtho (White Circle - Inside)", + Kokkos::RangePolicy(0, getBatchCount(start, end, offset)), + KOKKOS_LAMBDA(const int circle_task) { + const int i_r = start + circle_task * offset; + for (int i_theta = 0; i_theta < grid.ntheta(); i_theta++) { + nodeApplyAscOrthoCircleGiveInside(i_r, i_theta, grid, level_cache, DirBC_Interior, x, rhs, temp); + } + }); + Kokkos::fence(); } -} -} // namespace smoother_give + { + /* Outside White Section (Part 1)*/ + const int start = (grid.numberSmootherCircles() % 2 == 0) ? 1 : 0; + const int end = grid.numberSmootherCircles() + 1; + const int offset = 4; + Kokkos::parallel_for( + "SmootherGive: ApplyAscOrtho (White Circle - Outside: Part 1)", + Kokkos::RangePolicy(0, getBatchCount(start, end, offset)), + KOKKOS_LAMBDA(const int circle_task) { + const int i_r = start + circle_task * offset; + for (int i_theta = 0; i_theta < grid.ntheta(); i_theta++) { + nodeApplyAscOrthoCircleGiveOutside(i_r, i_theta, grid, level_cache, DirBC_Interior, x, rhs, temp); + } + }); + Kokkos::fence(); + } + + { + /* Outside White Section (Part 2)*/ + const int start = (grid.numberSmootherCircles() % 2 == 0) ? 3 : 2; + const int end = grid.numberSmootherCircles() + 1; + const int offset = 4; + Kokkos::parallel_for( + "SmootherGive: ApplyAscOrtho (White Circle - Outside: Part 2)", + Kokkos::RangePolicy(0, getBatchCount(start, end, offset)), + KOKKOS_LAMBDA(const int circle_task) { + const int i_r = start + circle_task * offset; + for (int i_theta = 0; i_theta < grid.ntheta(); i_theta++) { + nodeApplyAscOrthoCircleGiveOutside(i_r, i_theta, grid, level_cache, DirBC_Interior, x, rhs, temp); + } + }); + Kokkos::fence(); + } +} template -void SmootherGive::applyAscOrthoCircleSection(const int i_r, const SmootherColor smoother_color, - ConstVector x, ConstVector rhs, - Vector temp) +void SmootherGive::applyAscOrthoBlackRadialSection(ConstVector x, ConstVector rhs, + Vector temp) { - using smoother_give::nodeApplyAscOrthoCircleGive; + using smoother_give::nodeApplyAscOrthoRadialGiveInside; + using smoother_give::nodeApplyAscOrthoRadialGiveOutside; + + auto getBatchCount = [](int start, int end, int offset) { + if (start >= end) { + return 0; + } + return (end - start + offset - 1) / offset; + }; const PolarGrid& grid = Smoother::grid_; const LevelCacheType& level_cache = Smoother::level_cache_; const bool DirBC_Interior = Smoother::DirBC_Interior_; - assert(i_r >= 0 && i_r < grid.numberSmootherCircles() + 1); - - const double r = grid.radius(i_r); - - for (int i_theta = 0; i_theta < grid.ntheta(); i_theta++) { - const double theta = grid.theta(i_theta); - const int index = grid.index(i_r, i_theta); + /* ----------------------------------------------- */ + /* 3. Black-Radial update (u_br): */ + /* A_br * u_br = f_br − A_br^ortho * u_br^ortho */ + /* ----------------------------------------------- */ + + { + /* Inside Black Section */ + const int start = 0; + const int end = grid.ntheta(); + const int offset = 2; + + Kokkos::parallel_for( + "SmootherGive: ApplyAscOrtho (Black Radial - Inside)", + Kokkos::RangePolicy(0, getBatchCount(start, end, offset)), + KOKKOS_LAMBDA(const int radial_task) { + const int i_theta = start + radial_task * offset; + for (int i_r = grid.numberSmootherCircles() - 1; i_r < grid.nr(); i_r++) { + nodeApplyAscOrthoRadialGiveInside(i_r, i_theta, grid, level_cache, DirBC_Interior, x, rhs, temp); + } + }); + Kokkos::fence(); + } - double coeff_beta, arr, att, art, detDF; - level_cache.obtainValues(i_r, i_theta, index, r, theta, coeff_beta, arr, att, art, detDF); + { + /* Outside Black Section (Part 1) */ + const int start = 1; + const int end = grid.ntheta(); + const int offset = 4; + + Kokkos::parallel_for( + "SmootherGive: ApplyAscOrtho (Black Radial - Outside: Part 1)", + Kokkos::RangePolicy(0, getBatchCount(start, end, offset)), + KOKKOS_LAMBDA(const int radial_task) { + const int i_theta = start + radial_task * offset; + for (int i_r = grid.numberSmootherCircles(); i_r < grid.nr(); i_r++) { + nodeApplyAscOrthoRadialGiveOutside(i_r, i_theta, grid, level_cache, DirBC_Interior, x, rhs, temp); + } + }); + Kokkos::fence(); + } - // Apply Asc Ortho at the current node - nodeApplyAscOrthoCircleGive(i_r, i_theta, grid, DirBC_Interior, smoother_color, x, rhs, temp, arr, att, art, - detDF, coeff_beta); + { + /* Outside Black Section (Part 1) */ + const int start = 3; + const int end = grid.ntheta(); + const int offset = 4; + + Kokkos::parallel_for( + "SmootherGive: ApplyAscOrtho (Black Radial - Outside: Part 1)", + Kokkos::RangePolicy(0, getBatchCount(start, end, offset)), + KOKKOS_LAMBDA(const int radial_task) { + const int i_theta = start + radial_task * offset; + for (int i_r = grid.numberSmootherCircles(); i_r < grid.nr(); i_r++) { + nodeApplyAscOrthoRadialGiveOutside(i_r, i_theta, grid, level_cache, DirBC_Interior, x, rhs, temp); + } + }); + Kokkos::fence(); } } template -void SmootherGive::applyAscOrthoRadialSection(const int i_theta, const SmootherColor smoother_color, - ConstVector x, ConstVector rhs, - Vector temp) +void SmootherGive::applyAscOrthoWhiteRadialSection(ConstVector x, ConstVector rhs, + Vector temp) { - using smoother_give::nodeApplyAscOrthoRadialGive; + using smoother_give::nodeApplyAscOrthoRadialGiveInside; + using smoother_give::nodeApplyAscOrthoRadialGiveOutside; + + auto getBatchCount = [](int start, int end, int offset) { + if (start >= end) { + return 0; + } + return (end - start + offset - 1) / offset; + }; const PolarGrid& grid = Smoother::grid_; const LevelCacheType& level_cache = Smoother::level_cache_; const bool DirBC_Interior = Smoother::DirBC_Interior_; - const double theta = grid.theta(i_theta); - - /* We need to obtain left contributions from the circular section for AscOrtho. */ - /* !!! i_r = grid.numberSmootherCircles()-1 !!! */ - for (int i_r = grid.numberSmootherCircles() - 1; i_r < grid.nr(); i_r++) { - const double r = grid.radius(i_r); - const int index = grid.index(i_r, i_theta); + /* ----------------------------------------------- */ + /* 4. White-Radial update (u_wr): */ + /* A_wr * u_wr = f_wr − A_wr^ortho * u_wr^ortho */ + /* ----------------------------------------------- */ + + { + /* Inside White Section */ + const int start = 1; + const int end = grid.ntheta(); + const int offset = 2; + + Kokkos::parallel_for( + "SmootherGive: ApplyAscOrtho (White Radial - Inside)", + Kokkos::RangePolicy(0, getBatchCount(start, end, offset)), + KOKKOS_LAMBDA(const int radial_task) { + const int i_theta = start + radial_task * offset; + for (int i_r = grid.numberSmootherCircles() - 1; i_r < grid.nr(); i_r++) { + nodeApplyAscOrthoRadialGiveInside(i_r, i_theta, grid, level_cache, DirBC_Interior, x, rhs, temp); + } + }); + Kokkos::fence(); + } - double coeff_beta, arr, att, art, detDF; - level_cache.obtainValues(i_r, i_theta, index, r, theta, coeff_beta, arr, att, art, detDF); + { + /* Outside White Section (Part 1) */ + const int start = 0; + const int end = grid.ntheta(); + const int offset = 4; + + Kokkos::parallel_for( + "SmootherGive: ApplyAscOrtho (White Radial - Outside: Part 1)", + Kokkos::RangePolicy(0, getBatchCount(start, end, offset)), + KOKKOS_LAMBDA(const int radial_task) { + const int i_theta = start + radial_task * offset; + for (int i_r = grid.numberSmootherCircles(); i_r < grid.nr(); i_r++) { + nodeApplyAscOrthoRadialGiveOutside(i_r, i_theta, grid, level_cache, DirBC_Interior, x, rhs, temp); + } + }); + Kokkos::fence(); + } - // Apply Asc Ortho at the current node - nodeApplyAscOrthoRadialGive(i_r, i_theta, grid, DirBC_Interior, smoother_color, x, rhs, temp, arr, att, art, - detDF, coeff_beta); + { + /* Outside White Section (Part 1) */ + const int start = 2; + const int end = grid.ntheta(); + const int offset = 4; + + Kokkos::parallel_for( + "SmootherGive: ApplyAscOrtho (White Radial - Outside: Part 1)", + Kokkos::RangePolicy(0, getBatchCount(start, end, offset)), + KOKKOS_LAMBDA(const int radial_task) { + const int i_theta = start + radial_task * offset; + for (int i_r = grid.numberSmootherCircles(); i_r < grid.nr(); i_r++) { + nodeApplyAscOrthoRadialGiveOutside(i_r, i_theta, grid, level_cache, DirBC_Interior, x, rhs, temp); + } + }); + Kokkos::fence(); } -} \ No newline at end of file +} diff --git a/include/Smoother/SmootherGive/smootherGive.h b/include/Smoother/SmootherGive/smootherGive.h index 4265472b1..e3ac34689 100644 --- a/include/Smoother/SmootherGive/smootherGive.h +++ b/include/Smoother/SmootherGive/smootherGive.h @@ -130,6 +130,8 @@ class SmootherGive : public Smoother // Solver object (owns matrix if MUMPS, references if in-house solver). InnerBoundarySolver inner_boundary_solver_; + // Public is required as Cuda needs to be abe to get the address of functions enclosing lambda functions +public: /* -------------- */ /* Stencil access */ /* -------------- */ @@ -171,10 +173,10 @@ class SmootherGive : public Smoother // Compute temp = f_sc − A_sc^ortho * u_sc^ortho (precomputed right-hand side) // where x = u_sc and rhs = f_sc - void applyAscOrthoCircleSection(const int i_r, const SmootherColor smoother_color, ConstVector x, - ConstVector rhs, Vector temp); - void applyAscOrthoRadialSection(const int i_theta, const SmootherColor smoother_color, ConstVector x, - ConstVector rhs, Vector temp); + void applyAscOrthoBlackCircleSection(ConstVector x, ConstVector rhs, Vector temp); + void applyAscOrthoWhiteCircleSection(ConstVector x, ConstVector rhs, Vector temp); + void applyAscOrthoBlackRadialSection(ConstVector x, ConstVector rhs, Vector temp); + void applyAscOrthoWhiteRadialSection(ConstVector x, ConstVector rhs, Vector temp); /* ----------------- */ /* Line-wise solvers */ diff --git a/include/Smoother/SmootherGive/smootherGive.inl b/include/Smoother/SmootherGive/smootherGive.inl index 4081ef3b2..bd47d61f2 100644 --- a/include/Smoother/SmootherGive/smootherGive.inl +++ b/include/Smoother/SmootherGive/smootherGive.inl @@ -48,112 +48,15 @@ void SmootherGive::smoothing(Vector x, ConstVector::grid_; - const int num_omp_threads = Smoother::num_omp_threads_; - - /* Multi-threaded execution */ - const int num_smoother_circles = grid.numberSmootherCircles(); - const int num_radial_lines = grid.ntheta(); - - /* ----------------------------------------------- */ - /* 1. Black-Circle update (u_bc): */ - /* A_bc * u_bc = f_bc − A_bc^ortho * u_bc^ortho */ - /* ----------------------------------------------- */ -#pragma omp parallel num_threads(num_omp_threads) - { - /* Inside Black Section */ -#pragma omp for - for (int circle_task = 0; circle_task < num_smoother_circles; circle_task += 2) { - int i_r = num_smoother_circles - circle_task - 1; - applyAscOrthoCircleSection(i_r, SmootherColor::Black, x, rhs, temp); - } - /* Outside Black Section (Part 1)*/ -#pragma omp for - for (int circle_task = -1; circle_task < num_smoother_circles; circle_task += 4) { - int i_r = num_smoother_circles - circle_task - 1; - applyAscOrthoCircleSection(i_r, SmootherColor::Black, x, rhs, temp); - } - /* Outside Black Section (Part 2)*/ -#pragma omp for - for (int circle_task = 1; circle_task < num_smoother_circles; circle_task += 4) { - int i_r = num_smoother_circles - circle_task - 1; - applyAscOrthoCircleSection(i_r, SmootherColor::Black, x, rhs, temp); - } - } + applyAscOrthoBlackCircleSection(x, rhs, temp); solveBlackCircleSection(x, temp); - /* ----------------------------------------------- */ - /* 2. White-Circle update (u_wc): */ - /* A_wc * u_wc = f_wc − A_wc^ortho * u_wc^ortho */ - /* ----------------------------------------------- */ -#pragma omp parallel num_threads(num_omp_threads) - { - /* Inside White Section */ -#pragma omp for - for (int circle_task = 1; circle_task < num_smoother_circles; circle_task += 2) { - int i_r = num_smoother_circles - circle_task - 1; - applyAscOrthoCircleSection(i_r, SmootherColor::White, x, rhs, temp); - } - /* Outside White Section (Part 1)*/ -#pragma omp for - for (int circle_task = 0; circle_task < num_smoother_circles; circle_task += 4) { - int i_r = num_smoother_circles - circle_task - 1; - applyAscOrthoCircleSection(i_r, SmootherColor::White, x, rhs, temp); - } - /* Outside White Section (Part 2)*/ -#pragma omp for - for (int circle_task = 2; circle_task < num_smoother_circles; circle_task += 4) { - int i_r = num_smoother_circles - circle_task - 1; - applyAscOrthoCircleSection(i_r, SmootherColor::White, x, rhs, temp); - } - } + applyAscOrthoWhiteCircleSection(x, rhs, temp); solveWhiteCircleSection(x, temp); - /* ----------------------------------------------- */ - /* 3. Black-Radial update (u_br): */ - /* A_br * u_br = f_br − A_br^ortho * u_br^ortho */ - /* ----------------------------------------------- */ -#pragma omp parallel num_threads(num_omp_threads) - { - /* Inside Black Section */ -#pragma omp for - for (int i_theta = 0; i_theta < num_radial_lines; i_theta += 2) { - applyAscOrthoRadialSection(i_theta, SmootherColor::Black, x, rhs, temp); - } - /* Outside Black Section (Part 1) */ -#pragma omp for - for (int i_theta = 1; i_theta < num_radial_lines; i_theta += 4) { - applyAscOrthoRadialSection(i_theta, SmootherColor::Black, x, rhs, temp); - } - /* Outside Black Section (Part 2) */ -#pragma omp for - for (int i_theta = 3; i_theta < num_radial_lines; i_theta += 4) { - applyAscOrthoRadialSection(i_theta, SmootherColor::Black, x, rhs, temp); - } - } + applyAscOrthoBlackRadialSection(x, rhs, temp); solveBlackRadialSection(x, temp); - /* ----------------------------------------------- */ - /* 4. White-Radial update (u_wr): */ - /* A_wr * u_wr = f_wr − A_wr^ortho * u_wr^ortho */ - /* ----------------------------------------------- */ -#pragma omp parallel num_threads(num_omp_threads) - { - /* Inside Black Section */ -#pragma omp for - for (int i_theta = 1; i_theta < num_radial_lines; i_theta += 2) { - applyAscOrthoRadialSection(i_theta, SmootherColor::White, x, rhs, temp); - } - /* Outside Black Section (Part 1) */ -#pragma omp for - for (int i_theta = 0; i_theta < num_radial_lines; i_theta += 4) { - applyAscOrthoRadialSection(i_theta, SmootherColor::White, x, rhs, temp); - } - /* Outside Black Section (Part 2) */ -#pragma omp for - for (int i_theta = 2; i_theta < num_radial_lines; i_theta += 4) { - applyAscOrthoRadialSection(i_theta, SmootherColor::White, x, rhs, temp); - } - } + applyAscOrthoWhiteRadialSection(x, rhs, temp); solveWhiteRadialSection(x, temp); } From b04655deb5a15f887964eaf42e5cd1d1a78b6f3c Mon Sep 17 00:00:00 2001 From: julianlitz Date: Wed, 13 May 2026 19:48:17 +0200 Subject: [PATCH 02/10] hm --- .../Smoother/SmootherGive/applyAscOrtho.inl | 24 +++++++++---------- 1 file changed, 12 insertions(+), 12 deletions(-) diff --git a/include/Smoother/SmootherGive/applyAscOrtho.inl b/include/Smoother/SmootherGive/applyAscOrtho.inl index f220e4399..de38f82f4 100644 --- a/include/Smoother/SmootherGive/applyAscOrtho.inl +++ b/include/Smoother/SmootherGive/applyAscOrtho.inl @@ -94,10 +94,10 @@ static KOKKOS_INLINE_FUNCTION void nodeApplyAscOrthoCircleGiveInside(int i_r, in } template -static inline void nodeApplyAscOrthoCircleGiveOutside(int i_r, int i_theta, const PolarGrid& grid, - const LevelCacheType& level_cache, bool DirBC_Interior, - ConstVector& x, ConstVector& rhs, - Vector& result) +static KOKKOS_INLINE_FUNCTION void nodeApplyAscOrthoCircleGiveOutside(int i_r, int i_theta, const PolarGrid& grid, + const LevelCacheType& level_cache, + bool DirBC_Interior, ConstVector& x, + ConstVector& rhs, Vector& result) { assert(0 <= i_r && i_r <= grid.numberSmootherCircles()); @@ -146,10 +146,10 @@ static inline void nodeApplyAscOrthoCircleGiveOutside(int i_r, int i_theta, cons } template -static inline void nodeApplyAscOrthoRadialGiveInside(int i_r, int i_theta, const PolarGrid& grid, - const LevelCacheType& level_cache, bool DirBC_Interior, - ConstVector& x, ConstVector& rhs, - Vector& result) +static KOKKOS_INLINE_FUNCTION void nodeApplyAscOrthoRadialGiveInside(int i_r, int i_theta, const PolarGrid& grid, + const LevelCacheType& level_cache, + bool DirBC_Interior, ConstVector& x, + ConstVector& rhs, Vector& result) { assert(grid.numberSmootherCircles() - 1 <= i_r && i_r < grid.nr()); @@ -221,10 +221,10 @@ static inline void nodeApplyAscOrthoRadialGiveInside(int i_r, int i_theta, const } template -static inline void nodeApplyAscOrthoRadialGiveOutside(int i_r, int i_theta, const PolarGrid& grid, - const LevelCacheType& level_cache, bool DirBC_Interior, - ConstVector& x, ConstVector& rhs, - Vector& result) +static KOKKOS_INLINE_FUNCTION void nodeApplyAscOrthoRadialGiveOutside(int i_r, int i_theta, const PolarGrid& grid, + const LevelCacheType& level_cache, + bool DirBC_Interior, ConstVector& x, + ConstVector& rhs, Vector& result) { assert(grid.numberSmootherCircles() <= i_r && i_r < grid.nr()); From 865aa43a377affc162e9ee7174c30f4e0b907193 Mon Sep 17 00:00:00 2001 From: julianlitz Date: Wed, 13 May 2026 20:06:24 +0200 Subject: [PATCH 03/10] sus --- include/Smoother/SmootherGive/applyAscOrtho.inl | 6 +----- 1 file changed, 1 insertion(+), 5 deletions(-) diff --git a/include/Smoother/SmootherGive/applyAscOrtho.inl b/include/Smoother/SmootherGive/applyAscOrtho.inl index de38f82f4..8cce7b74f 100644 --- a/include/Smoother/SmootherGive/applyAscOrtho.inl +++ b/include/Smoother/SmootherGive/applyAscOrtho.inl @@ -32,8 +32,6 @@ static KOKKOS_INLINE_FUNCTION void nodeApplyAscOrthoCircleGiveInside(int i_r, in 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); @@ -119,8 +117,6 @@ static KOKKOS_INLINE_FUNCTION void nodeApplyAscOrthoCircleGiveOutside(int i_r, i 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); @@ -131,7 +127,7 @@ static KOKKOS_INLINE_FUNCTION void nodeApplyAscOrthoCircleGiveOutside(int i_r, i const int top = grid.index(i_r, i_theta_P1); /* Fill result(i-1,j) */ - if (i_r > 1 || !DirBC_Interior) { + if (i_r > 1 || (i_r == 1 && !DirBC_Interior)) { result[left] -= (-coeff1 * arr * x[center] /* Right */ - 0.25 * art * x[top] /* Top Right */ + 0.25 * art * x[bottom]); /* Bottom Right */ From f65cad7de082972418b539de5f702406099874d6 Mon Sep 17 00:00:00 2001 From: julianlitz Date: Wed, 13 May 2026 20:08:25 +0200 Subject: [PATCH 04/10] sus --- include/Smoother/SmootherGive/applyAscOrtho.inl | 17 ++++++++++------- 1 file changed, 10 insertions(+), 7 deletions(-) diff --git a/include/Smoother/SmootherGive/applyAscOrtho.inl b/include/Smoother/SmootherGive/applyAscOrtho.inl index 8cce7b74f..7ca02088b 100644 --- a/include/Smoother/SmootherGive/applyAscOrtho.inl +++ b/include/Smoother/SmootherGive/applyAscOrtho.inl @@ -110,30 +110,33 @@ static KOKKOS_INLINE_FUNCTION void nodeApplyAscOrthoCircleGiveOutside(int i_r, i level_cache.obtainValues(i_r, i_theta, center, radius, theta, coeff_beta, arr, att, art, detDF); if (0 <= i_r && 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 int i_theta_M1 = grid.wrapThetaIndex(i_theta - 1); const int i_theta_P1 = grid.wrapThetaIndex(i_theta + 1); - const int left = grid.index(i_r - 1, i_theta); - const int right = grid.index(i_r + 1, i_theta); const int bottom = grid.index(i_r, i_theta_M1); const int top = grid.index(i_r, i_theta_P1); /* Fill result(i-1,j) */ if (i_r > 1 || (i_r == 1 && !DirBC_Interior)) { + const double h1 = grid.radialSpacing(i_r - 1); + 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 */ } /* Fill result(i+1,j) */ if (i_r < grid.numberSmootherCircles() - 1) { + 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[right] -= (-coeff2 * arr * x[center] /* Left */ + 0.25 * art * x[top] /* Top Left */ - 0.25 * art * x[bottom]); /* Bottom Left */ From 97fa9d8d856303bc6a7026d05403217e570efa5e Mon Sep 17 00:00:00 2001 From: julianlitz Date: Wed, 13 May 2026 20:09:53 +0200 Subject: [PATCH 05/10] unused --- include/Smoother/SmootherGive/applyAscOrtho.inl | 2 -- 1 file changed, 2 deletions(-) diff --git a/include/Smoother/SmootherGive/applyAscOrtho.inl b/include/Smoother/SmootherGive/applyAscOrtho.inl index 7ca02088b..32832225c 100644 --- a/include/Smoother/SmootherGive/applyAscOrtho.inl +++ b/include/Smoother/SmootherGive/applyAscOrtho.inl @@ -243,8 +243,6 @@ static KOKKOS_INLINE_FUNCTION void nodeApplyAscOrthoRadialGiveOutside(int i_r, i double k1 = grid.angularSpacing(i_theta - 1); double k2 = grid.angularSpacing(i_theta); - double coeff1 = 0.5 * (k1 + k2) / h1; - double coeff2 = 0.5 * (k1 + k2) / h2; double coeff3 = 0.5 * (h1 + h2) / k1; double coeff4 = 0.5 * (h1 + h2) / k2; From aa09e0a8a3ba5949ca65263573f0be0cd4d22489 Mon Sep 17 00:00:00 2001 From: julianlitz Date: Wed, 13 May 2026 20:26:14 +0200 Subject: [PATCH 06/10] sus --- include/Smoother/SmootherGive/applyAscOrtho.inl | 6 +++--- 1 file changed, 3 insertions(+), 3 deletions(-) diff --git a/include/Smoother/SmootherGive/applyAscOrtho.inl b/include/Smoother/SmootherGive/applyAscOrtho.inl index 32832225c..d1e480d87 100644 --- a/include/Smoother/SmootherGive/applyAscOrtho.inl +++ b/include/Smoother/SmootherGive/applyAscOrtho.inl @@ -368,7 +368,7 @@ void SmootherGive::applyAscOrthoWhiteCircleSection(ConstVector::applyAscOrthoWhiteCircleSection(ConstVector::applyAscOrthoWhiteCircleSection(ConstVector Date: Wed, 13 May 2026 21:15:39 +0200 Subject: [PATCH 07/10] sanity --- .../Smoother/SmootherGive/applyAscOrtho.inl | 59 +++++++++++-------- 1 file changed, 33 insertions(+), 26 deletions(-) diff --git a/include/Smoother/SmootherGive/applyAscOrtho.inl b/include/Smoother/SmootherGive/applyAscOrtho.inl index d1e480d87..d94aba510 100644 --- a/include/Smoother/SmootherGive/applyAscOrtho.inl +++ b/include/Smoother/SmootherGive/applyAscOrtho.inl @@ -168,54 +168,61 @@ static KOKKOS_INLINE_FUNCTION void nodeApplyAscOrthoRadialGiveInside(int i_r, in const int bottom = grid.index(i_r, i_theta_M1); const int top = grid.index(i_r, i_theta_P1); - const int left = grid.index(i_r - 1, i_theta); - /* h1, k1, k2 are always valid at this point */ - const double h1 = grid.radialSpacing(i_r - 1); - 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 k1 = grid.angularSpacing(i_theta_M1); + const double k2 = grid.angularSpacing(i_theta); - /* ------------------- */ - /* Fill result(center) */ - /* ------------------- */ - if (grid.numberSmootherCircles() <= i_r && i_r < grid.nr() - 1) { - const double h2 = grid.radialSpacing(i_r); // safe: i_r < nr-1 + if (grid.numberSmootherCircles() <= i_r && i_r <= grid.nr() - 2) { + const double h1 = grid.radialSpacing(i_r); + const double h2 = grid.radialSpacing(i_r); const double coeff3 = 0.5 * (h1 + h2) / k1; const double coeff4 = 0.5 * (h1 + h2) / k2; result[center] -= (-coeff3 * att * x[bottom] /* Bottom */ - - coeff4 * att * x[top]); /* Top */ + - coeff4 * att * x[top]); /* Top */ } - /* ----------------- */ - /* Fill result(left) */ - /* ----------------- */ - if (grid.numberSmootherCircles() < i_r && i_r < grid.nr()) { - result[left] -= (-0.25 * art * x[top] /* Top Right */ + 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 */ } - /* ------------------- */ - /* Fill result(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 */ } - /* --------------- */ - /* Symmetry shifts */ - /* --------------- */ + if (i_r == grid.numberSmootherCircles() - 1) { + 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[right] -= (-coeff2 * arr * x[center]); /* Left */ + } + + if (i_r == grid.numberSmootherCircles()) { + const double h1 = grid.radialSpacing(i_r); + 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 */ + } + if (i_r == grid.nr() - 2) { 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 * rhs[right]; /* Right: Symmetry shift! */ + + const int right = grid.index(i_r + 1, i_theta); + result[center] -= (-coeff2 * arr * rhs[right]); /* Right: Symmetry shift! */ } if (i_r == grid.nr() - 1) { - result[left] -= -coeff1 * arr * rhs[center]; /* Right: Symmetry shift! */ + const double h1 = grid.radialSpacing(i_r); + const double coeff1 = 0.5 * (k1 + k2) / h1; + + const int left = grid.index(i_r - 1, i_theta); + result[left] -= (-coeff1 * arr * rhs[center]); /* Right: Symmetry shift! */ } } From 72c081c60549099a96c1ec4797247426faa0521f Mon Sep 17 00:00:00 2001 From: julianlitz Date: Wed, 13 May 2026 21:35:49 +0200 Subject: [PATCH 08/10] This is fine --- include/Smoother/SmootherGive/applyAscOrtho.inl | 6 +++--- 1 file changed, 3 insertions(+), 3 deletions(-) diff --git a/include/Smoother/SmootherGive/applyAscOrtho.inl b/include/Smoother/SmootherGive/applyAscOrtho.inl index d94aba510..2da54ec5f 100644 --- a/include/Smoother/SmootherGive/applyAscOrtho.inl +++ b/include/Smoother/SmootherGive/applyAscOrtho.inl @@ -173,7 +173,7 @@ static KOKKOS_INLINE_FUNCTION void nodeApplyAscOrthoRadialGiveInside(int i_r, in const double k2 = grid.angularSpacing(i_theta); if (grid.numberSmootherCircles() <= i_r && i_r <= grid.nr() - 2) { - const double h1 = grid.radialSpacing(i_r); + const double h1 = grid.radialSpacing(i_r - 1); const double h2 = grid.radialSpacing(i_r); const double coeff3 = 0.5 * (h1 + h2) / k1; const double coeff4 = 0.5 * (h1 + h2) / k2; @@ -202,7 +202,7 @@ static KOKKOS_INLINE_FUNCTION void nodeApplyAscOrthoRadialGiveInside(int i_r, in } if (i_r == grid.numberSmootherCircles()) { - const double h1 = grid.radialSpacing(i_r); + const double h1 = grid.radialSpacing(i_r - 1); const double coeff1 = 0.5 * (k1 + k2) / h1; const int left = grid.index(i_r - 1, i_theta); @@ -218,7 +218,7 @@ static KOKKOS_INLINE_FUNCTION void nodeApplyAscOrthoRadialGiveInside(int i_r, in } if (i_r == grid.nr() - 1) { - const double h1 = grid.radialSpacing(i_r); + const double h1 = grid.radialSpacing(i_r - 1); const double coeff1 = 0.5 * (k1 + k2) / h1; const int left = grid.index(i_r - 1, i_theta); From 1f02d97401a2c520d8dda0953aeddbfd2f09bee8 Mon Sep 17 00:00:00 2001 From: Julian Litz <91479202+julianlitz@users.noreply.github.com> Date: Fri, 15 May 2026 16:47:46 +0200 Subject: [PATCH 09/10] Update include/Smoother/SmootherGive/applyAscOrtho.inl Co-authored-by: Emily Bourne --- include/Smoother/SmootherGive/applyAscOrtho.inl | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/include/Smoother/SmootherGive/applyAscOrtho.inl b/include/Smoother/SmootherGive/applyAscOrtho.inl index 2da54ec5f..70bc59bbd 100644 --- a/include/Smoother/SmootherGive/applyAscOrtho.inl +++ b/include/Smoother/SmootherGive/applyAscOrtho.inl @@ -24,7 +24,7 @@ static KOKKOS_INLINE_FUNCTION void nodeApplyAscOrthoCircleGiveInside(int i_r, in /* -------------------- */ /* Node in the interior */ /* -------------------- */ - if (i_r > 0 && i_r < grid.numberSmootherCircles()) { + if (i_r > 0) { const double h1 = grid.radialSpacing(i_r - 1); const double h2 = grid.radialSpacing(i_r); const double k1 = grid.angularSpacing(i_theta - 1); From 574abad81987eea26419dcc64b8789bd3d318758 Mon Sep 17 00:00:00 2001 From: Julian Litz <91479202+julianlitz@users.noreply.github.com> Date: Fri, 15 May 2026 16:52:49 +0200 Subject: [PATCH 10/10] Update applyAscOrtho.inl --- .../Smoother/SmootherGive/applyAscOrtho.inl | 48 +++++++++---------- 1 file changed, 23 insertions(+), 25 deletions(-) diff --git a/include/Smoother/SmootherGive/applyAscOrtho.inl b/include/Smoother/SmootherGive/applyAscOrtho.inl index 70bc59bbd..f56c299c6 100644 --- a/include/Smoother/SmootherGive/applyAscOrtho.inl +++ b/include/Smoother/SmootherGive/applyAscOrtho.inl @@ -109,38 +109,36 @@ static KOKKOS_INLINE_FUNCTION void nodeApplyAscOrthoCircleGiveOutside(int i_r, i double coeff_beta, arr, att, art, detDF; level_cache.obtainValues(i_r, i_theta, center, radius, theta, coeff_beta, arr, att, art, detDF); - if (0 <= i_r && i_r <= grid.numberSmootherCircles()) { - const double k1 = grid.angularSpacing(i_theta - 1); - const double k2 = grid.angularSpacing(i_theta); + const double k1 = grid.angularSpacing(i_theta - 1); + 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 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 top = grid.index(i_r, i_theta_P1); + const int bottom = grid.index(i_r, i_theta_M1); + const int top = grid.index(i_r, i_theta_P1); - /* Fill result(i-1,j) */ - if (i_r > 1 || (i_r == 1 && !DirBC_Interior)) { - const double h1 = grid.radialSpacing(i_r - 1); - const double coeff1 = 0.5 * (k1 + k2) / h1; + /* Fill result(i-1,j) */ + if (i_r > 1 || (i_r == 1 && !DirBC_Interior)) { + const double h1 = grid.radialSpacing(i_r - 1); + const double coeff1 = 0.5 * (k1 + k2) / h1; - const int left = grid.index(i_r - 1, i_theta); + 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 */ - } - /* Fill result(i+1,j) */ - if (i_r < grid.numberSmootherCircles() - 1) { - const double h2 = grid.radialSpacing(i_r); - const double coeff2 = 0.5 * (k1 + k2) / h2; + result[left] -= (-coeff1 * arr * x[center] /* Right */ + - 0.25 * art * x[top] /* Top Right */ + + 0.25 * art * x[bottom]); /* Bottom Right */ + } + /* Fill result(i+1,j) */ + if (i_r < grid.numberSmootherCircles() - 1) { + 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); + 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 */ - } + result[right] -= (-coeff2 * arr * x[center] /* Left */ + + 0.25 * art * x[top] /* Top Left */ + - 0.25 * art * x[bottom]); /* Bottom Left */ } }