Skip to content

Add applySystemOperator(result, x) to Residual class for PCG support #186

New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

Sign up for GitHub

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Merged
julianlitz merged 1 commit into from
Feb 27, 2026
Merged

Add applySystemOperator(result, x) to Residual class for PCG support #186

Show file tree
Hide file tree
Changes from all commits
Commits
Show all changes
1 commit
Select commit
File filter

Filter by extension

Filter by extension
Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
2 changes: 2 additions & 0 deletions include/Residual/ResidualGive/residualGive.h
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -12,6 +12,8 @@ class ResidualGive : public Residual

void computeResidual(Vector<double> result, ConstVector<double> rhs, ConstVector<double> x) const override;

void applySystemOperator(Vector<double> result, ConstVector<double> x) const override;

private:
void applyCircleSection(const int i_r, Vector<double> result, ConstVector<double> x) const;
void applyRadialSection(const int i_theta, Vector<double> result, ConstVector<double> x) const;
Expand Down
7 changes: 4 additions & 3 deletions include/Residual/ResidualTake/residualTake.h
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -12,8 +12,9 @@ class ResidualTake : public Residual

void computeResidual(Vector<double> result, ConstVector<double> rhs, ConstVector<double> x) const override;

void applySystemOperator(Vector<double> result, ConstVector<double> x) const override;

private:
void applyCircleSection(const int i_r, Vector<double> result, ConstVector<double> rhs, ConstVector<double> x) const;
void applyRadialSection(const int i_theta, Vector<double> result, ConstVector<double> rhs,
ConstVector<double> x) const;
void applyCircleSection(const int i_r, Vector<double> result, ConstVector<double> x) const;
void applyRadialSection(const int i_theta, Vector<double> result, ConstVector<double> x) const;
};
2 changes: 2 additions & 0 deletions include/Residual/residual.h
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -25,6 +25,8 @@ class Residual

virtual void computeResidual(Vector<double> result, ConstVector<double> rhs, ConstVector<double> x) const = 0;

virtual void applySystemOperator(Vector<double> result, ConstVector<double> x) const = 0;

protected:
/* ------------------- */
/* Constructor members */
Expand Down
2 changes: 1 addition & 1 deletion src/CMakeLists.txt
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -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
)

Expand Down
46 changes: 23 additions & 23 deletions src/Residual/ResidualGive/applyAGive.cpp
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -34,30 +34,30 @@ 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 */
- coeff4 * att * x[top] /* Top */
+ ((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 */
Expand All @@ -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);
Expand All @@ -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 */
Expand Down Expand Up @@ -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 */
Expand All @@ -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 */
Expand Down Expand Up @@ -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 */
Expand All @@ -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 */
Expand Down Expand Up @@ -226,15 +226,15 @@ 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 */
- coeff4 * att * x[top] /* Top */
+ ((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 */
Expand All @@ -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 */
Expand All @@ -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);
Expand All @@ -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 */
Expand Down
24 changes: 20 additions & 4 deletions src/Residual/ResidualGive/residualGive.cpp
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -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<double> result, ConstVector<double> rhs, ConstVector<double> x) const
void ResidualGive::applySystemOperator(Vector<double> result, ConstVector<double> x) const
{
assert(result.size() == x.size());

Kokkos::deep_copy(result, rhs);
assign(result, 0.0);

/* Single-threaded execution */
if (num_omp_threads_ == 1) {
Expand Down Expand Up @@ -101,3 +101,19 @@ void ResidualGive::computeResidual(Vector<double> result, ConstVector<double> rh
}
}
// clang-format on

/* ------------------ */
/* result = rhs - A*x */
void ResidualGive::computeResidual(Vector<double> result, ConstVector<double> rhs, ConstVector<double> 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];
}
}
Loading