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

julianlitz · web-flow · commit 35f113ea5934 · 2026-02-27T16:47:19.000+01:00
diff --git a/include/Residual/ResidualGive/residualGive.h b/include/Residual/ResidualGive/residualGive.h
@@ -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;
diff --git a/include/Residual/ResidualTake/residualTake.h b/include/Residual/ResidualTake/residualTake.h
@@ -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;
 };
diff --git a/include/Residual/residual.h b/include/Residual/residual.h
@@ -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 */
diff --git 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
@@ -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 */
@@ -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,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 */
@@ -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
@@ -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) {
@@ -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];
+    }
+}
diff --git a/src/Residual/ResidualTake/applyATake.cpp b/src/Residual/ResidualTake/applyATake.cpp
diff --git a/src/Residual/ResidualTake/residualTake.cpp b/src/Residual/ResidualTake/residualTake.cpp

Original file line number	Diff line number	Diff line change
`@@ -104,7 +104,7 @@ set(RESIDUAL_SOURCES`
`104`	`104`	`${CMAKE_CURRENT_SOURCE_DIR}/Residual/ResidualGive/residualGive.cpp`
`105`	`105`
`106`	`106`	`# ResidualTake`
`107`		`- ${CMAKE_CURRENT_SOURCE_DIR}/Residual/ResidualTake/applyResidualTake.cpp`
	`107`	`+ ${CMAKE_CURRENT_SOURCE_DIR}/Residual/ResidualTake/applyATake.cpp`
`108`	`108`	`${CMAKE_CURRENT_SOURCE_DIR}/Residual/ResidualTake/residualTake.cpp`
`109`	`109`	`)`
`110`	`110`