From b3ef5237771cdcb4fb5a6357b1b08ff9f509299b Mon Sep 17 00:00:00 2001 From: julianlitz Date: Mon, 26 Jan 2026 23:30:52 +0100 Subject: [PATCH 1/6] Improve readability --- src/Residual/ResidualGive/applyAGive.cpp | 338 ++++++++++-------- src/Residual/ResidualGive/residualGive.cpp | 4 +- .../ResidualTake/applyResidualTake.cpp | 42 +-- src/Residual/ResidualTake/residualTake.cpp | 25 +- 4 files changed, 211 insertions(+), 198 deletions(-) diff --git a/src/Residual/ResidualGive/applyAGive.cpp b/src/Residual/ResidualGive/applyAGive.cpp index 3e2ed8e2..17a0cbf2 100644 --- a/src/Residual/ResidualGive/applyAGive.cpp +++ b/src/Residual/ResidualGive/applyAGive.cpp @@ -1,76 +1,91 @@ #include "../../../include/Residual/ResidualGive/residualGive.h" -#include "../../../include/common/geometry_helper.h" - inline void node_apply_a_give(int i_r, int i_theta, double r, double theta, const PolarGrid& grid, const LevelCache& level_cache, bool DirBC_Interior, Vector& result, ConstVector& x) { - const int global_index = grid.index(i_r, i_theta); + /* ---------------------------------------- */ + /* Compute or retrieve stencil coefficients */ + /* ---------------------------------------- */ + const int center = grid.index(i_r, i_theta); double coeff_beta, arr, att, art, detDF; - level_cache.obtainValues(i_r, i_theta, global_index, r, theta, coeff_beta, arr, att, art, detDF); + level_cache.obtainValues(i_r, i_theta, center, r, theta, coeff_beta, arr, att, art, detDF); + /* -------------------- */ /* Node in the interior */ /* -------------------- */ if (i_r > 1 && 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 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; + + 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 - 1); + const int top = grid.index(i_r, i_theta + 1); + /* Fill result(i,j) */ - result[grid.index(i_r, i_theta)] -= - (0.25 * (h1 + h2) * (k1 + k2) * coeff_beta * fabs(detDF) * x[grid.index(i_r, i_theta)] /* beta_{i,j} */ - - coeff1 * arr * x[grid.index(i_r - 1, i_theta)] /* Left */ - - coeff2 * arr * x[grid.index(i_r + 1, i_theta)] /* Right */ - - coeff3 * att * x[grid.index(i_r, i_theta - 1)] /* Bottom */ - - coeff4 * att * x[grid.index(i_r, i_theta + 1)] /* Top */ /* Center: (Left, Right, Bottom, Top) */ - + ((coeff1 + coeff2) * arr + (coeff3 + coeff4) * att) * x[grid.index(i_r, i_theta)]); + result[center] -= (0.25 * (h1 + h2) * (k1 + k2) * coeff_beta * fabs(detDF) * x[center] /* beta_{i,j} */ + - coeff1 * arr * x[left] /* Left */ + - coeff2 * arr * x[right] /* Right */ + - coeff3 * att * x[bottom] /* Bottom */ + - coeff4 * att * x[top] /* Top */ + + ((coeff1 + coeff2) * arr + (coeff3 + coeff4) * att) * + x[center]) /* Center: (Left, Right, Bottom, Top) */; /* Fill result(i-1,j) */ - result[grid.index(i_r - 1, i_theta)] -= (-coeff1 * arr * x[grid.index(i_r, i_theta)] /* Right */ - + coeff1 * arr * x[grid.index(i_r - 1, i_theta)] /* Center: (Right) */ - - 0.25 * art * x[grid.index(i_r, i_theta + 1)] /* Top Right */ - + 0.25 * art * x[grid.index(i_r, i_theta - 1)]); /* Bottom Right */ + 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 */ /* Fill result(i+1,j) */ - result[grid.index(i_r + 1, i_theta)] -= (-coeff2 * arr * x[grid.index(i_r, i_theta)] /* Left */ - + coeff2 * arr * x[grid.index(i_r + 1, i_theta)] /* Center: (Left) */ - + 0.25 * art * x[grid.index(i_r, i_theta + 1)] /* Top Left */ - - 0.25 * art * x[grid.index(i_r, i_theta - 1)]); /* Bottom Left */ + 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 */ /* Fill result(i,j-1) */ - result[grid.index(i_r, i_theta - 1)] -= (-coeff3 * att * x[grid.index(i_r, i_theta)] /* Top */ - + coeff3 * att * x[grid.index(i_r, i_theta - 1)] /* Center: (Top) */ - - 0.25 * art * x[grid.index(i_r + 1, i_theta)] /* Top Right */ - + 0.25 * art * x[grid.index(i_r - 1, i_theta)]); /* Top Left */ + 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 */ /* Fill result(i,j+1) */ - result[grid.index(i_r, i_theta + 1)] -= (-coeff4 * att * x[grid.index(i_r, i_theta)] /* Bottom */ - + coeff4 * att * x[grid.index(i_r, i_theta + 1)] /* Center: (Bottom) */ - + 0.25 * art * x[grid.index(i_r + 1, i_theta)] /* Bottom Right */ - - 0.25 * art * x[grid.index(i_r - 1, i_theta)]); /* Bottom Left */ - /* -------------------------- */ - /* Node on the inner boundary */ - /* -------------------------- */ + 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 */ } + /* -------------------------- */ + /* Node on the inner boundary */ + /* -------------------------- */ else if (i_r == 0) { /* ------------------------------------------------ */ /* Case 1: Dirichlet boundary on the inner boundary */ /* ------------------------------------------------ */ if (DirBC_Interior) { /* Fill result(i,j) */ - result[grid.index(i_r, i_theta)] -= x[grid.index(i_r, i_theta)]; + result[center] -= x[center]; + /* Give value to the interior nodes! */ - double h2 = grid.radialSpacing(i_r); - double k1 = grid.angularSpacing(i_theta - 1); - double k2 = grid.angularSpacing(i_theta); + 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; + + const int right = grid.index(i_r + 1, i_theta); + const int bottom = grid.index(i_r, i_theta - 1); + const int top = grid.index(i_r, i_theta + 1); + /* Fill result(i+1,j) */ - result[grid.index(i_r + 1, i_theta)] -= - (-coeff2 * arr * x[grid.index(i_r, i_theta)] /* Left */ - + coeff2 * arr * x[grid.index(i_r + 1, i_theta)] /* Center: (Left) */ - + 0.25 * art * x[grid.index(i_r, i_theta + 1)] /* Top Left */ - - 0.25 * art * x[grid.index(i_r, i_theta - 1)]); /* Bottom Left */ + 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 */ } else { /* ------------------------------------------------------------- */ @@ -79,156 +94,175 @@ inline void node_apply_a_give(int i_r, int i_theta, double r, double theta, cons /* 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 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 left = grid.index(i_r, i_theta + (grid.ntheta() / 2)); + const int right = grid.index(i_r + 1, i_theta); + const int bottom = grid.index(i_r, i_theta - 1); + const int top = grid.index(i_r, i_theta + 1); + /* Fill result(i,j) */ - result[grid.index(i_r, i_theta)] -= - (0.25 * (h1 + h2) * (k1 + k2) * coeff_beta * fabs(detDF) * x[grid.index(i_r, i_theta)] /* beta_{i,j} */ - - coeff1 * arr * x[grid.index(i_r, i_theta + (grid.ntheta() / 2))] /* Left */ - - coeff2 * arr * x[grid.index(i_r + 1, i_theta)] /* Right */ - - coeff3 * att * x[grid.index(i_r, i_theta - 1)] /* Bottom */ - - coeff4 * att * x[grid.index(i_r, i_theta + 1)] /* Top */ - + ((coeff1 + coeff2) * arr + (coeff3 + coeff4) * att) * - x[grid.index(i_r, i_theta)]); /* Center: (Left, Right, Bottom, Top) */ + result[center] -= (0.25 * (h1 + h2) * (k1 + k2) * coeff_beta * fabs(detDF) * x[center] /* beta_{i,j} */ + - coeff1 * arr * x[left] /* Left */ + - coeff2 * arr * x[right] /* Right */ + - coeff3 * att * x[bottom] /* Bottom */ + - coeff4 * att * x[top] /* Top */ + + ((coeff1 + coeff2) * arr + (coeff3 + coeff4) * att) * + x[center]); /* Center: (Left, Right, Bottom, Top) */ /* Fill result(i-1,j) */ /* From view the view of the across origin node, the directions are roatated by 180 degrees in the stencil! */ - result[grid.index(i_r, i_theta + (grid.ntheta() / 2))] -= - (-coeff1 * arr * x[grid.index(i_r, i_theta)] /* Right -> Left */ - + coeff1 * arr * - x[grid.index(i_r, i_theta + (grid.ntheta() / 2))]); /* Center: (Right) -> Center: (Left)*/ - /* + 0.25 * art * x[grid.index(i_r,i_theta+1)]; // Top Right -> Bottom Left: REMOVED DUE TO ARTIFICAL 7 POINT STENCIL */ - /* - 0.25 * art * x[grid.index(i_r,i_theta-1)]; // Bottom Right -> Top Left: REMOVED DUE TO ARTIFICAL 7 POINT 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 */ /* Fill result(i+1,j) */ - result[grid.index(i_r + 1, i_theta)] -= - (-coeff2 * arr * x[grid.index(i_r, i_theta)] /* Left */ - + coeff2 * arr * x[grid.index(i_r + 1, i_theta)] /* Center: (Left) */ - + 0.25 * art * x[grid.index(i_r, i_theta + 1)] /* Top Left */ - - 0.25 * art * x[grid.index(i_r, i_theta - 1)]); /* Bottom Left */ + 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 */ /* Fill result(i,j-1) */ - result[grid.index(i_r, i_theta - 1)] -= - (-coeff3 * att * x[grid.index(i_r, i_theta)] /* Top */ - + coeff3 * att * x[grid.index(i_r, i_theta - 1)] /* Center: (Top) */ - - 0.25 * art * x[grid.index(i_r + 1, i_theta)]); /* Top Right */ - /* + 0.25 * art * x[grid.index(i_r-1,i_theta)]; // Top Left: REMOVED DUE TO ARTIFICAL 7 POINT STENCIL */ + 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 */ /* Fill result(i,j+1) */ - result[grid.index(i_r, i_theta + 1)] -= - (-coeff4 * att * x[grid.index(i_r, i_theta)] /* Bottom */ - + coeff4 * att * x[grid.index(i_r, i_theta + 1)] /* Center: (Bottom) */ - + 0.25 * art * x[grid.index(i_r + 1, i_theta)]); /* Bottom Right */ - /* - 0.25 * art * x[grid.index(i_r-1,i_theta)]; // Bottom Left: REMOVED DUE TO ARTIFICAL 7 POINT STENCIL */ + 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 */ } - /* ------------------------------- */ - /* Node next to the inner boundary */ - /* ------------------------------- */ } + /* ------------------------------- */ + /* Node next to the inner boundary */ + /* ------------------------------- */ else if (i_r == 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); + 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; + + 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 - 1); + const int top = grid.index(i_r, i_theta + 1); + /* Fill result(i,j) */ - result[grid.index(i_r, i_theta)] -= - (0.25 * (h1 + h2) * (k1 + k2) * coeff_beta * fabs(detDF) * x[grid.index(i_r, i_theta)] /* beta_{i,j} */ - - coeff1 * arr * x[grid.index(i_r - 1, i_theta)] /* Left */ - - coeff2 * arr * x[grid.index(i_r + 1, i_theta)] /* Right */ - - coeff3 * att * x[grid.index(i_r, i_theta - 1)] /* Bottom */ - - coeff4 * att * x[grid.index(i_r, i_theta + 1)] /* Top */ - + ((coeff1 + coeff2) * arr + (coeff3 + coeff4) * att) * - x[grid.index(i_r, i_theta)]); /* Center: (Left, Right, Bottom, Top) */ + result[center] -= (0.25 * (h1 + h2) * (k1 + k2) * coeff_beta * fabs(detDF) * x[center] /* beta_{i,j} */ + - coeff1 * arr * x[left] /* Left */ + - coeff2 * arr * x[right] /* Right */ + - coeff3 * att * x[bottom] /* Bottom */ + - coeff4 * att * x[top] /* Top */ + + ((coeff1 + coeff2) * arr + (coeff3 + coeff4) * att) * + x[center]); /* Center: (Left, Right, Bottom, Top) */ /* Fill result(i-1,j) */ if (!DirBC_Interior) { /* Don't give to the inner dirichlet boundary! */ - result[grid.index(i_r - 1, i_theta)] -= - (-coeff1 * arr * x[grid.index(i_r, i_theta)] /* Right */ - + coeff1 * arr * x[grid.index(i_r - 1, i_theta)] /* Center: (Right) */ - - 0.25 * art * x[grid.index(i_r, i_theta + 1)] /* Top Right */ - + 0.25 * art * x[grid.index(i_r, i_theta - 1)]); /* Bottom Right */ + 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 */ } /* Fill result(i+1,j) */ - result[grid.index(i_r + 1, i_theta)] -= (-coeff2 * arr * x[grid.index(i_r, i_theta)] /* Left */ - + coeff2 * arr * x[grid.index(i_r + 1, i_theta)] /* Center: (Left) */ - + 0.25 * art * x[grid.index(i_r, i_theta + 1)] /* Top Left */ - - 0.25 * art * x[grid.index(i_r, i_theta - 1)]); /* Bottom Left */ + 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 */ /* Fill result(i,j-1) */ - result[grid.index(i_r, i_theta - 1)] -= (-coeff3 * att * x[grid.index(i_r, i_theta)] /* Top */ - + coeff3 * att * x[grid.index(i_r, i_theta - 1)] /* Center: (Top) */ - - 0.25 * art * x[grid.index(i_r + 1, i_theta)] /* Top Right */ - + 0.25 * art * x[grid.index(i_r - 1, i_theta)]); /* Top Left */ + 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 */ /* Fill result(i,j+1) */ - result[grid.index(i_r, i_theta + 1)] -= (-coeff4 * att * x[grid.index(i_r, i_theta)] /* Bottom */ - + coeff4 * att * x[grid.index(i_r, i_theta + 1)] /* Center: (Bottom) */ - + 0.25 * art * x[grid.index(i_r + 1, i_theta)] /* Bottom Right */ - - 0.25 * art * x[grid.index(i_r - 1, i_theta)]); /* Bottom Left */ - /* ------------------------------- */ - /* Node next to the outer boundary */ - /* ------------------------------- */ + 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 */ } + /* ------------------------------- */ + /* Node next to the outer 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 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; + + 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 - 1); + const int top = grid.index(i_r, i_theta + 1); + /* Fill result(i,j) */ - result[grid.index(i_r, i_theta)] -= - (0.25 * (h1 + h2) * (k1 + k2) * coeff_beta * fabs(detDF) * x[grid.index(i_r, i_theta)] /* beta_{i,j} */ - - coeff1 * arr * x[grid.index(i_r - 1, i_theta)] /* Left */ - - coeff2 * arr * x[grid.index(i_r + 1, i_theta)] /* Right */ - - coeff3 * att * x[grid.index(i_r, i_theta - 1)] /* Bottom */ - - coeff4 * att * x[grid.index(i_r, i_theta + 1)] /* Top */ - + ((coeff1 + coeff2) * arr + (coeff3 + coeff4) * att) * - x[grid.index(i_r, i_theta)]); /* Center: (Left, Right, Bottom, Top) */ + result[center] -= (0.25 * (h1 + h2) * (k1 + k2) * coeff_beta * fabs(detDF) * x[center] /* beta_{i,j} */ + - coeff1 * arr * x[left] /* Left */ + - coeff2 * arr * x[right] /* Right */ + - coeff3 * att * x[bottom] /* Bottom */ + - coeff4 * att * x[top] /* Top */ + + ((coeff1 + coeff2) * arr + (coeff3 + coeff4) * att) * + x[center]); /* Center: (Left, Right, Bottom, Top) */ /* Fill result(i-1,j) */ - result[grid.index(i_r - 1, i_theta)] -= (-coeff1 * arr * x[grid.index(i_r, i_theta)] /* Right */ - + coeff1 * arr * x[grid.index(i_r - 1, i_theta)] /* Center: (Right) */ - - 0.25 * art * x[grid.index(i_r, i_theta + 1)] /* Top Right */ - + 0.25 * art * x[grid.index(i_r, i_theta - 1)]); /* Bottom Right */ + 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 */ /* Don't give to the outer dirichlet boundary! */ /* Fill result(i+1,j) */ - /* result[grid.index(i_r+1,i_theta)] -= ( */ - /* - coeff2 * arr * x[grid.index(i_r,i_theta)] // Left */ - /* + coeff2 * arr * x[grid.index(i_r+1,i_theta)] // Center: (Left) */ - /* + 0.25 * art * x[grid.index(i_r,i_theta+1)] // Top Left */ - /* - 0.25 * art * x[grid.index(i_r,i_theta-1)] ); // Bottom Left */ + /* result[grid.index(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 */ /* Fill result(i,j-1) */ - result[grid.index(i_r, i_theta - 1)] -= (-coeff3 * att * x[grid.index(i_r, i_theta)] /* Top */ - + coeff3 * att * x[grid.index(i_r, i_theta - 1)] /* Center: (Top) */ - - 0.25 * art * x[grid.index(i_r + 1, i_theta)] /* Top Right */ - + 0.25 * art * x[grid.index(i_r - 1, i_theta)]); /* Top Left */ + 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 */ /* Fill result(i,j+1) */ - result[grid.index(i_r, i_theta + 1)] -= (-coeff4 * att * x[grid.index(i_r, i_theta)] /* Bottom */ - + coeff4 * att * x[grid.index(i_r, i_theta + 1)] /* Center: (Bottom) */ - + 0.25 * art * x[grid.index(i_r + 1, i_theta)] /* Bottom Right */ - - 0.25 * art * x[grid.index(i_r - 1, i_theta)]); /* Bottom Left */ - /* ----------------------------- */ - /* Node on to the outer boundary */ - /* ----------------------------- */ + 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 */ } + /* ----------------------------- */ + /* Node on to the outer boundary */ + /* ----------------------------- */ else if (i_r == grid.nr() - 1) { /* Fill result of (i,j) */ - result[grid.index(i_r, i_theta)] -= x[grid.index(i_r, i_theta)]; + result[center] -= x[center]; + /* Give value to the interior nodes! */ - double h1 = grid.radialSpacing(i_r - 1); - double k1 = grid.angularSpacing(i_theta - 1); - double k2 = grid.angularSpacing(i_theta); + 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; + + const int left = grid.index(i_r - 1, i_theta); + const int bottom = grid.index(i_r, i_theta - 1); + const int top = grid.index(i_r, i_theta + 1); + /* Fill result(i-1,j) */ - result[grid.index(i_r - 1, i_theta)] -= (-coeff1 * arr * x[grid.index(i_r, i_theta)] /* Right */ - + coeff1 * arr * x[grid.index(i_r - 1, i_theta)] /* Center: (Right) */ - - 0.25 * art * x[grid.index(i_r, i_theta + 1)] /* Top Right */ - + 0.25 * art * x[grid.index(i_r, i_theta - 1)]); /* Bottom Right */ + 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 */ } } @@ -237,7 +271,6 @@ void ResidualGive::applyCircleSection(const int i_r, Vector result, Cons 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); - node_apply_a_give(i_r, i_theta, r, theta, grid_, level_cache_, DirBC_Interior_, result, x); } } @@ -247,7 +280,6 @@ void ResidualGive::applyRadialSection(const int i_theta, Vector result, const double theta = grid_.theta(i_theta); for (int i_r = grid_.numberSmootherCircles(); i_r < grid_.nr(); i_r++) { const double r = grid_.radius(i_r); - node_apply_a_give(i_r, i_theta, r, theta, grid_, level_cache_, DirBC_Interior_, result, x); } } diff --git a/src/Residual/ResidualGive/residualGive.cpp b/src/Residual/ResidualGive/residualGive.cpp index 9e6dc047..eb308998 100644 --- a/src/Residual/ResidualGive/residualGive.cpp +++ b/src/Residual/ResidualGive/residualGive.cpp @@ -18,8 +18,8 @@ void ResidualGive::computeResidual(Vector result, ConstVector rh Kokkos::deep_copy(result, rhs); + /* Single-threaded execution */ if (num_omp_threads_ == 1) { - /* Single-threaded execution */ for (int i_r = 0; i_r < grid_.numberSmootherCircles(); i_r++) { applyCircleSection(i_r, result, x); } @@ -27,8 +27,8 @@ void ResidualGive::computeResidual(Vector result, ConstVector rh applyRadialSection(i_theta, result, x); } } + /* Multi-threaded execution */ else { - /* Multi-threaded execution */ const int num_circle_tasks = grid_.numberSmootherCircles(); const int additional_radial_tasks = grid_.ntheta() % 3; const int num_radial_tasks = grid_.ntheta() - additional_radial_tasks; diff --git a/src/Residual/ResidualTake/applyResidualTake.cpp b/src/Residual/ResidualTake/applyResidualTake.cpp index b4443dad..e202f91d 100644 --- a/src/Residual/ResidualTake/applyResidualTake.cpp +++ b/src/Residual/ResidualTake/applyResidualTake.cpp @@ -19,18 +19,15 @@ inline void node_apply_residual_take(int i_r, int i_theta, const PolarGrid& grid 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 bottom_left = grid.index(i_r - 1, i_theta_M1); + const int bottom_left = grid.index(i_r - 1, i_theta - 1); 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_left = grid.index(i_r - 1, i_theta + 1); + const int bottom = grid.index(i_r, i_theta - 1); 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 top = grid.index(i_r, i_theta + 1); + const int bottom_right = grid.index(i_r + 1, i_theta - 1); const int right = grid.index(i_r + 1, i_theta); - const int top_right = grid.index(i_r + 1, i_theta_P1); + const int top_right = grid.index(i_r + 1, i_theta + 1); result[center] = rhs[center] - @@ -46,10 +43,10 @@ inline void node_apply_residual_take(int i_r, int i_theta, const PolarGrid& grid + 0.25 * (art[left] + art[top]) * x[top_left] /* Top Left */ - 0.25 * (art[right] + art[top]) * x[top_right] /* Top Right */ ); - /* -------------------------- */ - /* Node on the inner boundary */ - /* -------------------------- */ } + /* -------------------------- */ + /* Node on the inner boundary */ + /* -------------------------- */ else if (i_r == 0) { /* ------------------------------------------------ */ /* Case 1: Dirichlet boundary on the inner boundary */ @@ -75,17 +72,13 @@ inline void node_apply_residual_take(int i_r, int i_theta, const PolarGrid& grid 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 int bottom = grid.index(i_r, i_theta_M1); + const int left = grid.index(i_r, i_theta + (grid.ntheta() / 2)); + const int bottom = grid.index(i_r, i_theta - 1); 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 top = grid.index(i_r, i_theta + 1); + const int bottom_right = grid.index(i_r + 1, i_theta - 1); const int right = grid.index(i_r + 1, i_theta); - const int top_right = grid.index(i_r + 1, i_theta_P1); + const int top_right = grid.index(i_r + 1, i_theta + 1); result[center] = rhs[center] - @@ -98,15 +91,14 @@ inline void node_apply_residual_take(int i_r, int i_theta, const PolarGrid& grid /* - 0.25 * (art[left] + art[bottom]) * x[bottom_left] // Bottom Left: REMOVED DUE TO ARTIFICAL 7 POINT STENCIL */ + 0.25 * (art[right] + art[bottom]) * x[bottom_right] /* Bottom Right */ - /* + 0.25 * (art[left] + art[top]) * x[top_left] // Top Left: REMOVED DUE TO ARTIFICAL 7 POINT STENCIL */ - 0.25 * (art[right] + art[top]) * x[top_right] /* Top Right */ ); } - /* ----------------------------- */ - /* Node on to the outer boundary */ - /* ----------------------------- */ } + /* ----------------------------- */ + /* Node on to the outer boundary */ + /* ----------------------------- */ else if (i_r == grid.nr() - 1) { /* Fill result of (i,j) */ const int center = grid.index(i_r, i_theta); diff --git a/src/Residual/ResidualTake/residualTake.cpp b/src/Residual/ResidualTake/residualTake.cpp index 3cbbd13e..ca10c140 100644 --- a/src/Residual/ResidualTake/residualTake.cpp +++ b/src/Residual/ResidualTake/residualTake.cpp @@ -18,30 +18,19 @@ void ResidualTake::computeResidual(Vector result, ConstVector rh assert(level_cache_.cacheDensityProfileCoefficients()); assert(level_cache_.cacheDomainGeometry()); - if (num_omp_threads_ == 1) { - /* Single-threaded execution */ + #pragma omp parallel num_threads(num_omp_threads_) + { + /* Circle Section */ + #pragma omp for nowait for (int i_r = 0; i_r < grid_.numberSmootherCircles(); i_r++) { applyCircleSection(i_r, result, rhs, x); } + + /* Radial Section */ + #pragma omp for nowait for (int i_theta = 0; i_theta < grid_.ntheta(); i_theta++) { applyRadialSection(i_theta, result, rhs, x); } } - else { - /* Multi-threaded execution */ - #pragma omp parallel num_threads(num_omp_threads_) - { - /* Circle Section */ - #pragma omp for nowait - for (int i_r = 0; i_r < grid_.numberSmootherCircles(); i_r++) { - applyCircleSection(i_r, result, rhs, x); - } - /* Radial Section */ - #pragma omp for nowait - for (int i_theta = 0; i_theta < grid_.ntheta(); i_theta++) { - applyRadialSection(i_theta, result, rhs, x); - } - } - } } // clang-format on From 131b278ba249ef1b077a53f3452d651bbdc4e489 Mon Sep 17 00:00:00 2001 From: julianlitz Date: Mon, 26 Jan 2026 23:39:05 +0100 Subject: [PATCH 2/6] happy little accident --- src/Residual/ResidualGive/applyAGive.cpp | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/Residual/ResidualGive/applyAGive.cpp b/src/Residual/ResidualGive/applyAGive.cpp index 17a0cbf2..afc5c8cd 100644 --- a/src/Residual/ResidualGive/applyAGive.cpp +++ b/src/Residual/ResidualGive/applyAGive.cpp @@ -224,7 +224,7 @@ inline void node_apply_a_give(int i_r, int i_theta, double r, double theta, cons + 0.25 * art * x[bottom]); /* Bottom Right */ /* Don't give to the outer dirichlet boundary! */ /* Fill result(i+1,j) */ - /* result[grid.index(right] -= ( */ + /* result[right] -= ( */ /* - coeff2 * arr * x[center] // Left */ /* + coeff2 * arr * x[right] // Center: (Left) */ /* + 0.25 * art * x[top] // Top Left */ From 7552a1a60f848bb81867e490548958d71bb89010 Mon Sep 17 00:00:00 2001 From: julianlitz Date: Mon, 26 Jan 2026 23:46:55 +0100 Subject: [PATCH 3/6] Use std::fabs --- src/Residual/ResidualGive/applyAGive.cpp | 8 ++++---- 1 file changed, 4 insertions(+), 4 deletions(-) diff --git a/src/Residual/ResidualGive/applyAGive.cpp b/src/Residual/ResidualGive/applyAGive.cpp index afc5c8cd..17dd666d 100644 --- a/src/Residual/ResidualGive/applyAGive.cpp +++ b/src/Residual/ResidualGive/applyAGive.cpp @@ -31,7 +31,7 @@ inline void node_apply_a_give(int i_r, int i_theta, double r, double theta, cons const int top = grid.index(i_r, i_theta + 1); /* Fill result(i,j) */ - result[center] -= (0.25 * (h1 + h2) * (k1 + k2) * coeff_beta * fabs(detDF) * x[center] /* beta_{i,j} */ + result[center] -= (0.25 * (h1 + h2) * (k1 + k2) * coeff_beta * std::fabs(detDF) * x[center] /* beta_{i,j} */ - coeff1 * arr * x[left] /* Left */ - coeff2 * arr * x[right] /* Right */ - coeff3 * att * x[bottom] /* Bottom */ @@ -110,7 +110,7 @@ inline void node_apply_a_give(int i_r, int i_theta, double r, double theta, cons const int top = grid.index(i_r, i_theta + 1); /* Fill result(i,j) */ - result[center] -= (0.25 * (h1 + h2) * (k1 + k2) * coeff_beta * fabs(detDF) * x[center] /* beta_{i,j} */ + result[center] -= (0.25 * (h1 + h2) * (k1 + k2) * coeff_beta * std::fabs(detDF) * x[center] /* beta_{i,j} */ - coeff1 * arr * x[left] /* Left */ - coeff2 * arr * x[right] /* Right */ - coeff3 * att * x[bottom] /* Bottom */ @@ -160,7 +160,7 @@ inline void node_apply_a_give(int i_r, int i_theta, double r, double theta, cons const int top = grid.index(i_r, i_theta + 1); /* Fill result(i,j) */ - result[center] -= (0.25 * (h1 + h2) * (k1 + k2) * coeff_beta * fabs(detDF) * x[center] /* beta_{i,j} */ + result[center] -= (0.25 * (h1 + h2) * (k1 + k2) * coeff_beta * std::fabs(detDF) * x[center] /* beta_{i,j} */ - coeff1 * arr * x[left] /* Left */ - coeff2 * arr * x[right] /* Right */ - coeff3 * att * x[bottom] /* Bottom */ @@ -210,7 +210,7 @@ inline void node_apply_a_give(int i_r, int i_theta, double r, double theta, cons const int top = grid.index(i_r, i_theta + 1); /* Fill result(i,j) */ - result[center] -= (0.25 * (h1 + h2) * (k1 + k2) * coeff_beta * fabs(detDF) * x[center] /* beta_{i,j} */ + result[center] -= (0.25 * (h1 + h2) * (k1 + k2) * coeff_beta * std::fabs(detDF) * x[center] /* beta_{i,j} */ - coeff1 * arr * x[left] /* Left */ - coeff2 * arr * x[right] /* Right */ - coeff3 * att * x[bottom] /* Bottom */ From cbc4182e78268e942e5530bd39845192d5243ae8 Mon Sep 17 00:00:00 2001 From: Julian Litz <91479202+julianlitz@users.noreply.github.com> Date: Mon, 26 Jan 2026 23:50:39 +0100 Subject: [PATCH 4/6] Update applyResidualTake.cpp --- src/Residual/ResidualTake/applyResidualTake.cpp | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/src/Residual/ResidualTake/applyResidualTake.cpp b/src/Residual/ResidualTake/applyResidualTake.cpp index e202f91d..181bec45 100644 --- a/src/Residual/ResidualTake/applyResidualTake.cpp +++ b/src/Residual/ResidualTake/applyResidualTake.cpp @@ -31,7 +31,7 @@ inline void node_apply_residual_take(int i_r, int i_theta, const PolarGrid& grid result[center] = rhs[center] - - (0.25 * (h1 + h2) * (k1 + k2) * coeff_beta[center] * fabs(detDF[center]) * x[center] /* beta_{i,j} */ + (0.25 * (h1 + h2) * (k1 + k2) * coeff_beta[center] * std::fabs(detDF[center]) * x[center] /* beta_{i,j} */ - coeff1 * (arr[center] + arr[left]) * (x[left] - x[center]) /* Left - Center: (Left) */ - coeff2 * (arr[center] + arr[right]) * (x[right] - x[center]) /* Right - Center: (Right) */ @@ -82,7 +82,7 @@ inline void node_apply_residual_take(int i_r, int i_theta, const PolarGrid& grid result[center] = rhs[center] - - (0.25 * (h1 + h2) * (k1 + k2) * coeff_beta[center] * fabs(detDF[center]) * x[center] /* beta_{i,j} */ + (0.25 * (h1 + h2) * (k1 + k2) * coeff_beta[center] * std::fabs(detDF[center]) * x[center] /* beta_{i,j} */ - coeff1 * (arr[center] + arr[left]) * (x[left] - x[center]) /* Left - Center: (Left) */ - coeff2 * (arr[center] + arr[right]) * (x[right] - x[center]) /* Right - Center: (Right) */ From 362a2ece77b278cde5009bd7bc2807f38fc00ae7 Mon Sep 17 00:00:00 2001 From: Julian Litz <91479202+julianlitz@users.noreply.github.com> Date: Tue, 27 Jan 2026 00:03:40 +0100 Subject: [PATCH 5/6] Update applyResidualTake.cpp --- src/Residual/ResidualTake/applyResidualTake.cpp | 3 ++- 1 file changed, 2 insertions(+), 1 deletion(-) diff --git a/src/Residual/ResidualTake/applyResidualTake.cpp b/src/Residual/ResidualTake/applyResidualTake.cpp index 181bec45..880b16b0 100644 --- a/src/Residual/ResidualTake/applyResidualTake.cpp +++ b/src/Residual/ResidualTake/applyResidualTake.cpp @@ -82,7 +82,8 @@ inline void node_apply_residual_take(int i_r, int i_theta, const PolarGrid& grid result[center] = rhs[center] - - (0.25 * (h1 + h2) * (k1 + k2) * coeff_beta[center] * std::fabs(detDF[center]) * x[center] /* beta_{i,j} */ + (0.25 * (h1 + h2) * (k1 + k2) * coeff_beta[center] * std::fabs(detDF[center]) * + x[center] /* beta_{i,j} */ - coeff1 * (arr[center] + arr[left]) * (x[left] - x[center]) /* Left - Center: (Left) */ - coeff2 * (arr[center] + arr[right]) * (x[right] - x[center]) /* Right - Center: (Right) */ From c69f02e52242f7933b6109a5318069c7fed88e03 Mon Sep 17 00:00:00 2001 From: julianlitz Date: Tue, 27 Jan 2026 13:19:09 +0100 Subject: [PATCH 6/6] Explicit wrapping --- src/Residual/ResidualGive/applyAGive.cpp | 45 +++++++++++++------ .../ResidualTake/applyResidualTake.cpp | 29 +++++++----- 2 files changed, 50 insertions(+), 24 deletions(-) diff --git a/src/Residual/ResidualGive/applyAGive.cpp b/src/Residual/ResidualGive/applyAGive.cpp index 17dd666d..b250bbc5 100644 --- a/src/Residual/ResidualGive/applyAGive.cpp +++ b/src/Residual/ResidualGive/applyAGive.cpp @@ -25,10 +25,13 @@ inline void node_apply_a_give(int i_r, int i_theta, double r, double theta, cons 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 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 - 1); - const int top = grid.index(i_r, i_theta + 1); + 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] -= (0.25 * (h1 + h2) * (k1 + k2) * coeff_beta * std::fabs(detDF) * x[center] /* beta_{i,j} */ @@ -77,9 +80,12 @@ inline void node_apply_a_give(int i_r, int i_theta, double r, double theta, cons 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 right = grid.index(i_r + 1, i_theta); - const int bottom = grid.index(i_r, i_theta - 1); - const int top = grid.index(i_r, i_theta + 1); + 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) */ result[right] -= (-coeff2 * arr * x[center] /* Left */ @@ -104,10 +110,14 @@ inline void node_apply_a_give(int i_r, int i_theta, double r, double theta, cons double coeff3 = 0.5 * (h1 + h2) / k1; double coeff4 = 0.5 * (h1 + h2) / k2; - const int left = grid.index(i_r, i_theta + (grid.ntheta() / 2)); + const int i_theta_M1 = grid.wrapThetaIndex(i_theta - 1); + const int i_theta_P1 = grid.wrapThetaIndex(i_theta + 1); + const int i_theta_Across = grid.wrapThetaIndex(i_theta + grid.ntheta() / 2); + + const int left = grid.index(i_r, i_theta_Across); const int right = grid.index(i_r + 1, i_theta); - const int bottom = grid.index(i_r, i_theta - 1); - const int top = grid.index(i_r, i_theta + 1); + 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] -= (0.25 * (h1 + h2) * (k1 + k2) * coeff_beta * std::fabs(detDF) * x[center] /* beta_{i,j} */ @@ -154,10 +164,13 @@ inline void node_apply_a_give(int i_r, int i_theta, double r, double theta, cons 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 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 - 1); - const int top = grid.index(i_r, i_theta + 1); + 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] -= (0.25 * (h1 + h2) * (k1 + k2) * coeff_beta * std::fabs(detDF) * x[center] /* beta_{i,j} */ @@ -204,10 +217,13 @@ inline void node_apply_a_give(int i_r, int i_theta, double r, double theta, cons 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 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 - 1); - const int top = grid.index(i_r, i_theta + 1); + 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] -= (0.25 * (h1 + h2) * (k1 + k2) * coeff_beta * std::fabs(detDF) * x[center] /* beta_{i,j} */ @@ -254,9 +270,12 @@ inline void node_apply_a_give(int i_r, int i_theta, double r, double theta, cons double coeff1 = 0.5 * (k1 + k2) / h1; + const int i_theta_M1 = grid.wrapThetaIndex(i_theta - 1); + const int i_theta_P1 = grid.wrapThetaIndex(i_theta + 1); + const int left = grid.index(i_r - 1, i_theta); - const int bottom = grid.index(i_r, i_theta - 1); - const int top = grid.index(i_r, i_theta + 1); + 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) */ result[left] -= (-coeff1 * arr * x[center] /* Right */ diff --git a/src/Residual/ResidualTake/applyResidualTake.cpp b/src/Residual/ResidualTake/applyResidualTake.cpp index 880b16b0..11140ff8 100644 --- a/src/Residual/ResidualTake/applyResidualTake.cpp +++ b/src/Residual/ResidualTake/applyResidualTake.cpp @@ -19,15 +19,18 @@ inline void node_apply_residual_take(int i_r, int i_theta, const PolarGrid& grid double coeff3 = 0.5 * (h1 + h2) / k1; double coeff4 = 0.5 * (h1 + h2) / k2; - const int bottom_left = grid.index(i_r - 1, i_theta - 1); + const int i_theta_M1 = grid.wrapThetaIndex(i_theta - 1); + const int i_theta_P1 = grid.wrapThetaIndex(i_theta + 1); + + const int bottom_left = grid.index(i_r - 1, i_theta_M1); const int left = grid.index(i_r - 1, i_theta); - const int top_left = grid.index(i_r - 1, i_theta + 1); - const int bottom = grid.index(i_r, i_theta - 1); + 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 + 1); - const int bottom_right = grid.index(i_r + 1, i_theta - 1); + 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 + 1); + const int top_right = grid.index(i_r + 1, i_theta_P1); result[center] = rhs[center] - @@ -72,13 +75,17 @@ inline void node_apply_residual_take(int i_r, int i_theta, const PolarGrid& grid double coeff3 = 0.5 * (h1 + h2) / k1; double coeff4 = 0.5 * (h1 + h2) / k2; - const int left = grid.index(i_r, i_theta + (grid.ntheta() / 2)); - const int bottom = grid.index(i_r, 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 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 + 1); - const int bottom_right = grid.index(i_r + 1, i_theta - 1); + 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 + 1); + const int top_right = grid.index(i_r + 1, i_theta_P1); result[center] = rhs[center] -