diff --git a/src/Residual/ResidualGive/applyAGive.cpp b/src/Residual/ResidualGive/applyAGive.cpp index 3e2ed8e2..b250bbc5 100644 --- a/src/Residual/ResidualGive/applyAGive.cpp +++ b/src/Residual/ResidualGive/applyAGive.cpp @@ -1,76 +1,97 @@ #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 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,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 * std::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 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_M1); + const int top = grid.index(i_r, i_theta_P1); + /* 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 +100,188 @@ 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 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_M1); + const int top = grid.index(i_r, i_theta_P1); + /* 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 * std::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 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,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 * std::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 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,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 * std::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[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 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_M1); + const int top = grid.index(i_r, i_theta_P1); + /* 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 +290,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 +299,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..11140ff8 100644 --- a/src/Residual/ResidualTake/applyResidualTake.cpp +++ b/src/Residual/ResidualTake/applyResidualTake.cpp @@ -34,7 +34,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) */ @@ -46,10 +46,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 */ @@ -89,7 +89,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] * 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) */ @@ -98,15 +99,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