diff --git a/include/Residual/ResidualGive/residualGive.h b/include/Residual/ResidualGive/residualGive.h index a26e5f8f..d489cc2d 100644 --- a/include/Residual/ResidualGive/residualGive.h +++ b/include/Residual/ResidualGive/residualGive.h @@ -12,6 +12,8 @@ class ResidualGive : public Residual void computeResidual(Vector result, ConstVector rhs, ConstVector x) const override; + void applySystemOperator(Vector result, ConstVector x) const override; + private: void applyCircleSection(const int i_r, Vector result, ConstVector x) const; void applyRadialSection(const int i_theta, Vector result, ConstVector x) const; diff --git a/include/Residual/ResidualTake/residualTake.h b/include/Residual/ResidualTake/residualTake.h index 5c3a3648..4907e2dc 100644 --- a/include/Residual/ResidualTake/residualTake.h +++ b/include/Residual/ResidualTake/residualTake.h @@ -12,8 +12,9 @@ class ResidualTake : public Residual void computeResidual(Vector result, ConstVector rhs, ConstVector x) const override; + void applySystemOperator(Vector result, ConstVector x) const override; + private: - void applyCircleSection(const int i_r, Vector result, ConstVector rhs, ConstVector x) const; - void applyRadialSection(const int i_theta, Vector result, ConstVector rhs, - ConstVector x) const; + void applyCircleSection(const int i_r, Vector result, ConstVector x) const; + void applyRadialSection(const int i_theta, Vector result, ConstVector x) const; }; \ No newline at end of file diff --git a/include/Residual/residual.h b/include/Residual/residual.h index ba85e9bf..6317f302 100644 --- a/include/Residual/residual.h +++ b/include/Residual/residual.h @@ -25,6 +25,8 @@ class Residual virtual void computeResidual(Vector result, ConstVector rhs, ConstVector x) const = 0; + virtual void applySystemOperator(Vector result, ConstVector x) const = 0; + protected: /* ------------------- */ /* Constructor members */ diff --git a/src/CMakeLists.txt b/src/CMakeLists.txt index 0b00e56f..1427f30b 100644 --- a/src/CMakeLists.txt +++ b/src/CMakeLists.txt @@ -104,7 +104,7 @@ set(RESIDUAL_SOURCES ${CMAKE_CURRENT_SOURCE_DIR}/Residual/ResidualGive/residualGive.cpp # ResidualTake - ${CMAKE_CURRENT_SOURCE_DIR}/Residual/ResidualTake/applyResidualTake.cpp + ${CMAKE_CURRENT_SOURCE_DIR}/Residual/ResidualTake/applyATake.cpp ${CMAKE_CURRENT_SOURCE_DIR}/Residual/ResidualTake/residualTake.cpp ) diff --git a/src/Residual/ResidualGive/applyAGive.cpp b/src/Residual/ResidualGive/applyAGive.cpp index 15193fa4..61b83e69 100644 --- a/src/Residual/ResidualGive/applyAGive.cpp +++ b/src/Residual/ResidualGive/applyAGive.cpp @@ -34,7 +34,7 @@ static inline void node_apply_a_give(int i_r, int i_theta, double r, double thet 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} */ + 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 */ @@ -42,22 +42,22 @@ static inline void node_apply_a_give(int i_r, int i_theta, double r, double thet + ((coeff1 + coeff2) * arr + (coeff3 + coeff4) * att) * x[center]) /* Center: (Left, Right, Bottom, Top) */; /* Fill result(i-1,j) */ - result[left] -= (-coeff1 * arr * x[center] /* 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[right] -= (-coeff2 * arr * x[center] /* 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[bottom] -= (-coeff3 * att * x[center] /* Top */ + 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[top] -= (-coeff4 * att * x[center] /* Bottom */ + 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 */ @@ -71,7 +71,7 @@ static inline void node_apply_a_give(int i_r, int i_theta, double r, double thet /* ------------------------------------------------ */ if (DirBC_Interior) { /* Fill result(i,j) */ - result[center] -= x[center]; + result[center] += x[center]; /* Give value to the interior nodes! */ double h2 = grid.radialSpacing(i_r); @@ -88,7 +88,7 @@ static inline void node_apply_a_give(int i_r, int i_theta, double r, double thet const int top = grid.index(i_r, i_theta_P1); /* Fill result(i+1,j) */ - result[right] -= (-coeff2 * arr * x[center] /* 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 */ @@ -120,7 +120,7 @@ static inline void node_apply_a_give(int i_r, int i_theta, double r, double thet 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} */ + 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 */ @@ -129,22 +129,22 @@ static inline void node_apply_a_give(int i_r, int i_theta, double r, double thet 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[left] -= (-coeff1 * arr * x[center] /* Right -> Left */ + 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[right] -= (-coeff2 * arr * x[center] /* 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[bottom] -= (-coeff3 * att * x[center] /* Top */ + 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[top] -= (-coeff4 * att * x[center] /* Bottom */ + 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 */ @@ -173,7 +173,7 @@ static inline void node_apply_a_give(int i_r, int i_theta, double r, double thet 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} */ + 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 */ @@ -182,23 +182,23 @@ static inline void node_apply_a_give(int i_r, int i_theta, double r, double thet x[center]); /* Center: (Left, Right, Bottom, Top) */ /* Fill result(i-1,j) */ if (!DirBC_Interior) { /* Don't give to the inner dirichlet boundary! */ - result[left] -= (-coeff1 * arr * x[center] /* 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[right] -= (-coeff2 * arr * x[center] /* 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[bottom] -= (-coeff3 * att * x[center] /* Top */ + 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[top] -= (-coeff4 * att * x[center] /* Bottom */ + 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 */ @@ -226,7 +226,7 @@ static inline void node_apply_a_give(int i_r, int i_theta, double r, double thet 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} */ + 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 */ @@ -234,7 +234,7 @@ static inline void node_apply_a_give(int i_r, int i_theta, double r, double thet + ((coeff1 + coeff2) * arr + (coeff3 + coeff4) * att) * x[center]); /* Center: (Left, Right, Bottom, Top) */ /* Fill result(i-1,j) */ - result[left] -= (-coeff1 * arr * x[center] /* 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 */ @@ -246,12 +246,12 @@ static inline void node_apply_a_give(int i_r, int i_theta, double r, double thet /* + 0.25 * art * x[top] // Top Left */ /* - 0.25 * art * x[bottom]); // Bottom Left */ /* Fill result(i,j-1) */ - result[bottom] -= (-coeff3 * att * x[center] /* Top */ + 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[top] -= (-coeff4 * att * x[center] /* Bottom */ + 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 */ @@ -261,7 +261,7 @@ static inline void node_apply_a_give(int i_r, int i_theta, double r, double thet /* ----------------------------- */ else if (i_r == grid.nr() - 1) { /* Fill result of (i,j) */ - result[center] -= x[center]; + result[center] += x[center]; /* Give value to the interior nodes! */ double h1 = grid.radialSpacing(i_r - 1); @@ -278,7 +278,7 @@ static inline void node_apply_a_give(int i_r, int i_theta, double r, double thet const int top = grid.index(i_r, i_theta_P1); /* Fill result(i-1,j) */ - result[left] -= (-coeff1 * arr * x[center] /* 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 */ diff --git a/src/Residual/ResidualGive/residualGive.cpp b/src/Residual/ResidualGive/residualGive.cpp index eb308998..8ae4c38f 100644 --- a/src/Residual/ResidualGive/residualGive.cpp +++ b/src/Residual/ResidualGive/residualGive.cpp @@ -8,15 +8,15 @@ ResidualGive::ResidualGive(const PolarGrid& grid, const LevelCache& level_cache, { } -/* ------------------ */ -/* result = rhs - A*x */ +/* ------------ */ +/* result = A*x */ // clang-format off -void ResidualGive::computeResidual(Vector result, ConstVector rhs, ConstVector x) const +void ResidualGive::applySystemOperator(Vector result, ConstVector x) const { assert(result.size() == x.size()); - Kokkos::deep_copy(result, rhs); + assign(result, 0.0); /* Single-threaded execution */ if (num_omp_threads_ == 1) { @@ -101,3 +101,19 @@ void ResidualGive::computeResidual(Vector result, ConstVector rh } } // clang-format on + +/* ------------------ */ +/* result = rhs - A*x */ +void ResidualGive::computeResidual(Vector result, ConstVector rhs, ConstVector x) const +{ + assert(result.size() == x.size()); + + applySystemOperator(result, x); + + // Subtract A*x from rhs to get the residual. + const int n = result.size(); +#pragma omp parallel for num_threads(num_omp_threads_) + for (int i = 0; i < n; i++) { + result[i] = rhs[i] - result[i]; + } +} \ No newline at end of file diff --git a/src/Residual/ResidualTake/applyResidualTake.cpp b/src/Residual/ResidualTake/applyATake.cpp similarity index 61% rename from src/Residual/ResidualTake/applyResidualTake.cpp rename to src/Residual/ResidualTake/applyATake.cpp index 9c4a2417..ba62de85 100644 --- a/src/Residual/ResidualTake/applyResidualTake.cpp +++ b/src/Residual/ResidualTake/applyATake.cpp @@ -1,10 +1,9 @@ #include "../../../include/Residual/ResidualTake/residualTake.h" -static inline void node_apply_residual_take(int i_r, int i_theta, const PolarGrid& grid, bool DirBC_Interior, - Vector& result, ConstVector& rhs, ConstVector& x, - ConstVector& arr, ConstVector& att, - ConstVector& art, ConstVector& detDF, - ConstVector& coeff_beta) +static inline void node_apply_a_take(int i_r, int i_theta, const PolarGrid& grid, bool DirBC_Interior, + Vector& result, ConstVector& x, ConstVector& arr, + ConstVector& att, ConstVector& art, ConstVector& detDF, + ConstVector& coeff_beta) { /* -------------------- */ /* Node in the interior */ @@ -34,19 +33,17 @@ static inline void node_apply_residual_take(int i_r, int i_theta, const PolarGri const int top_right = grid.index(i_r + 1, i_theta_P1); result[center] = - rhs[center] - - (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) */ - - coeff3 * (att[center] + att[bottom]) * (x[bottom] - x[center]) /* Bottom - Center: (Bottom) */ - - coeff4 * (att[center] + att[top]) * (x[top] - x[center]) /* Top - Center: (Top) */ - - - 0.25 * (art[left] + art[bottom]) * x[bottom_left] /* Bottom Left */ - + 0.25 * (art[right] + art[bottom]) * x[bottom_right] /* Bottom Right */ - + 0.25 * (art[left] + art[top]) * x[top_left] /* Top Left */ - - 0.25 * (art[right] + art[top]) * x[top_right] /* Top Right */ - ); + 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) */ + - coeff3 * (att[center] + att[bottom]) * (x[bottom] - x[center]) /* Bottom - Center: (Bottom) */ + - coeff4 * (att[center] + att[top]) * (x[top] - x[center]) /* Top - Center: (Top) */ + + - 0.25 * (art[left] + art[bottom]) * x[bottom_left] /* Bottom Left */ + + 0.25 * (art[right] + art[bottom]) * x[bottom_right] /* Bottom Right */ + + 0.25 * (art[left] + art[top]) * x[top_left] /* Top Left */ + - 0.25 * (art[right] + art[top]) * x[top_right]; /* Top Right */ } /* -------------------------- */ /* Node on the inner boundary */ @@ -57,7 +54,7 @@ static inline void node_apply_residual_take(int i_r, int i_theta, const PolarGri /* ------------------------------------------------ */ if (DirBC_Interior) { const int center = grid.index(i_r, i_theta); - result[center] = rhs[center] - x[center]; + result[center] = x[center]; } else { /* ------------------------------------------------------------- */ @@ -89,20 +86,18 @@ static inline void node_apply_residual_take(int i_r, int i_theta, const PolarGri const int top_right = grid.index(i_r + 1, i_theta_P1); result[center] = - rhs[center] - - (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) */ - - coeff3 * (att[center] + att[bottom]) * (x[bottom] - x[center]) /* Bottom - Center: (Bottom) */ - - coeff4 * (att[center] + att[top]) * (x[top] - x[center]) /* Top - Center: (Top) */ - - /* - 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 */ - ); + 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) */ + - coeff3 * (att[center] + att[bottom]) * (x[bottom] - x[center]) /* Bottom - Center: (Bottom) */ + - coeff4 * (att[center] + att[top]) * (x[top] - x[center]) /* Top - Center: (Top) */ + + /* - 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 */ } } /* ----------------------------- */ @@ -111,12 +106,11 @@ static inline void node_apply_residual_take(int i_r, int i_theta, const PolarGri else if (i_r == grid.nr() - 1) { /* Fill result of (i,j) */ const int center = grid.index(i_r, i_theta); - result[center] = rhs[center] - x[center]; + result[center] = x[center]; } } -void ResidualTake::applyCircleSection(const int i_r, Vector result, ConstVector rhs, - ConstVector x) const +void ResidualTake::applyCircleSection(const int i_r, Vector result, ConstVector x) const { assert(level_cache_.cacheDensityProfileCoefficients()); assert(level_cache_.cacheDomainGeometry()); @@ -128,13 +122,11 @@ void ResidualTake::applyCircleSection(const int i_r, Vector result, Cons ConstVector coeff_beta = level_cache_.coeff_beta(); for (int i_theta = 0; i_theta < grid_.ntheta(); i_theta++) { - node_apply_residual_take(i_r, i_theta, grid_, DirBC_Interior_, result, rhs, x, arr, att, art, detDF, - coeff_beta); + node_apply_a_take(i_r, i_theta, grid_, DirBC_Interior_, result, x, arr, att, art, detDF, coeff_beta); } } -void ResidualTake::applyRadialSection(const int i_theta, Vector result, ConstVector rhs, - ConstVector x) const +void ResidualTake::applyRadialSection(const int i_theta, Vector result, ConstVector x) const { assert(level_cache_.cacheDensityProfileCoefficients()); assert(level_cache_.cacheDomainGeometry()); @@ -146,7 +138,6 @@ void ResidualTake::applyRadialSection(const int i_theta, Vector result, ConstVector coeff_beta = level_cache_.coeff_beta(); for (int i_r = grid_.numberSmootherCircles(); i_r < grid_.nr(); i_r++) { - node_apply_residual_take(i_r, i_theta, grid_, DirBC_Interior_, result, rhs, x, arr, att, art, detDF, - coeff_beta); + node_apply_a_take(i_r, i_theta, grid_, DirBC_Interior_, result, x, arr, att, art, detDF, coeff_beta); } } diff --git a/src/Residual/ResidualTake/residualTake.cpp b/src/Residual/ResidualTake/residualTake.cpp index ca10c140..c75df393 100644 --- a/src/Residual/ResidualTake/residualTake.cpp +++ b/src/Residual/ResidualTake/residualTake.cpp @@ -7,30 +7,43 @@ ResidualTake::ResidualTake(const PolarGrid& grid, const LevelCache& level_cache, { } -/* ------------------ */ -/* result = rhs - A*x */ - -// clang-format off -void ResidualTake::computeResidual(Vector result, ConstVector rhs, ConstVector x) const +/* ------------ */ +/* result = A*x */ +void ResidualTake::applySystemOperator(Vector result, ConstVector x) const { assert(result.size() == x.size()); assert(level_cache_.cacheDensityProfileCoefficients()); assert(level_cache_.cacheDomainGeometry()); - #pragma omp parallel num_threads(num_omp_threads_) +#pragma omp parallel num_threads(num_omp_threads_) { - /* Circle Section */ - #pragma omp for nowait +/* Circle Section */ +#pragma omp for nowait for (int i_r = 0; i_r < grid_.numberSmootherCircles(); i_r++) { - applyCircleSection(i_r, result, rhs, x); + applyCircleSection(i_r, result, x); } - - /* Radial Section */ - #pragma omp for nowait + +/* Radial Section */ +#pragma omp for nowait for (int i_theta = 0; i_theta < grid_.ntheta(); i_theta++) { - applyRadialSection(i_theta, result, rhs, x); + applyRadialSection(i_theta, result, x); } } } -// clang-format on + +/* ------------------ */ +/* result = rhs - A*x */ +void ResidualTake::computeResidual(Vector result, ConstVector rhs, ConstVector x) const +{ + assert(result.size() == x.size()); + + applySystemOperator(result, x); + + // Subtract A*x from rhs to get the residual. + const int n = result.size(); +#pragma omp parallel for num_threads(num_omp_threads_) + for (int i = 0; i < n; i++) { + result[i] = rhs[i] - result[i]; + } +} \ No newline at end of file