Substitute COO_Mumps_Solver in SmootheGive

include/Smoother/SmootherGive/applyAscOrtho.inl

@@ -434,6 +434,7 @@ void SmootherGive<DomainGeometry>::applyAscOrthoCircleSection(const int i_r, con
 
     const PolarGrid& grid                         = Smoother<DomainGeometry>::grid_;
     const LevelCache<DomainGeometry>& level_cache = Smoother<DomainGeometry>::level_cache_;
+    const bool DirBC_Interior                     = Smoother<DomainGeometry>::DirBC_Interior_;
 
     assert(i_r >= 0 && i_r < grid.numberSmootherCircles() + 1);
 
@@ -447,8 +448,8 @@ void SmootherGive<DomainGeometry>::applyAscOrthoCircleSection(const int i_r, con
         level_cache.obtainValues(i_r, i_theta, index, r, theta, coeff_beta, arr, att, art, detDF);
 
         // Apply Asc Ortho at the current node
-        nodeApplyAscOrthoCircleGive(i_r, i_theta, grid, Smoother<DomainGeometry>::DirBC_Interior_, smoother_color, x,
-                                    rhs, temp, arr, att, art, detDF, coeff_beta);
+        nodeApplyAscOrthoCircleGive(i_r, i_theta, grid, DirBC_Interior, smoother_color, x, rhs, temp, arr, att, art,
+                                    detDF, coeff_beta);
     }
 }
 
@@ -461,6 +462,7 @@ void SmootherGive<DomainGeometry>::applyAscOrthoRadialSection(const int i_theta,
 
     const PolarGrid& grid                         = Smoother<DomainGeometry>::grid_;
     const LevelCache<DomainGeometry>& level_cache = Smoother<DomainGeometry>::level_cache_;
+    const bool DirBC_Interior                     = Smoother<DomainGeometry>::DirBC_Interior_;
 
     const double theta = grid.theta(i_theta);
 
@@ -474,7 +476,7 @@ void SmootherGive<DomainGeometry>::applyAscOrthoRadialSection(const int i_theta,
         level_cache.obtainValues(i_r, i_theta, index, r, theta, coeff_beta, arr, att, art, detDF);
 
         // Apply Asc Ortho at the current node
-        nodeApplyAscOrthoRadialGive(i_r, i_theta, grid, Smoother<DomainGeometry>::DirBC_Interior_, smoother_color, x,
-                                    rhs, temp, arr, att, art, detDF, coeff_beta);
+        nodeApplyAscOrthoRadialGive(i_r, i_theta, grid, DirBC_Interior, smoother_color, x, rhs, temp, arr, att, art,
+                                    detDF, coeff_beta);
     }
 }

include/Smoother/SmootherGive/buildInnerBoundaryAsc.inl

@@ -0,0 +1,290 @@
+#pragma once
+
+namespace smoother_give
+{
+
+#ifdef GMGPOLAR_USE_MUMPS
+// When using the MUMPS solver, the matrix is assembled in COO format.
+static inline void update_CSR_COO_MatrixElement(SparseMatrixCOO<double>& matrix, int ptr, int offset, int row,
+                                                int column, double value)
+{
+    matrix.row_index(ptr + offset) = row;
+    matrix.col_index(ptr + offset) = column;
+    matrix.value(ptr + offset) += value;
+}
+#else
+// When using the in-house solver, the matrix is stored in CSR format.
+static inline void update_CSR_COO_MatrixElement(SparseMatrixCSR<double>& matrix, int ptr, int offset, int row,
+                                                int column, double value)
+{
+    matrix.row_nz_index(row, offset) = column;
+    matrix.row_nz_entry(row, offset) += value;
+}
+#endif
+
+} // namespace smoother_give
+
+template <concepts::DomainGeometry DomainGeometry>
+void SmootherGive<DomainGeometry>::nodeBuildInteriorBoundarySolverMatrix_i_r_0(int i_theta, const PolarGrid& grid,
+                                                                               bool DirBC_Interior,
+                                                                               InnerBoundaryMatrix& matrix, double arr,
+                                                                               double att, double art, double detDF,
+                                                                               double coeff_beta)
+{
+    using smoother_give::update_CSR_COO_MatrixElement;
+
+    assert(i_theta >= 0 && i_theta < grid.ntheta());
+
+    int ptr, offset;
+    int row, column;
+    double value;
+
+    const int i_r = 0;
+
+    /* ------------------------------------------------ */
+    /* Case 1: Dirichlet boundary on the inner boundary */
+    /* ------------------------------------------------ */
+    if (DirBC_Interior) {
+        /* Fill result(i,j) */
+        const int center_index = i_theta;
+
+        /* Fill matrix row of (i,j) */
+        row = center_index;
+        ptr = getCircleAscIndex(i_r, i_theta);
+
+        const Stencil& CenterStencil = getStencil(i_r);
+
+        offset = CenterStencil[StencilPosition::Center];
+        column = center_index;
+        value  = 1.0;
+        update_CSR_COO_MatrixElement(matrix, ptr, offset, row, column, value);
+    }
+    else {
+        /* ------------------------------------------------------------- */
+        /* Case 2: Across origin discretization on the interior boundary */
+        /* ------------------------------------------------------------- */
+        // h1 gets replaced with 2 * R0.
+        // (i_r-1,i_theta) gets replaced with (i_r, i_theta + (grid.ntheta()>>1)).
+        // Some more adjustments from the changing the 9-point stencil to the artifical 7-point stencil.
+        const double h1 = 2 * grid.radius(0);
+        const double h2 = grid.radialSpacing(i_r);
+        const double k1 = grid.angularSpacing(i_theta - 1);
+        const double k2 = grid.angularSpacing(i_theta);
+
+        const double coeff1 = 0.5 * (k1 + k2) / h1;
+        const double coeff2 = 0.5 * (k1 + k2) / h2;
+        const double coeff3 = 0.5 * (h1 + h2) / k1;
+        const double coeff4 = 0.5 * (h1 + h2) / k2;
+
+        /* -| x | o | x | */
+        /* -|   |   |   | */
+        /* -| O | o | o | */
+        /* -|   |   |   | */
+        /* -| x | o | x | */
+        const int i_theta_M1           = grid.wrapThetaIndex(i_theta - 1);
+        const int i_theta_P1           = grid.wrapThetaIndex(i_theta + 1);
+        const int i_theta_AcrossOrigin = grid.wrapThetaIndex(i_theta + grid.ntheta() / 2);
+
+        const int center_index = i_theta;
+        const int left_index   = i_theta_AcrossOrigin;
+        const int bottom_index = i_theta_M1;
+        const int top_index    = i_theta_P1;
+
+        const int center_nz_index = getCircleAscIndex(i_r, i_theta);
+        const int bottom_nz_index = getCircleAscIndex(i_r, i_theta_M1);
+        const int top_nz_index    = getCircleAscIndex(i_r, i_theta_P1);
+        const int left_nz_index   = getCircleAscIndex(i_r, i_theta_AcrossOrigin);
+
+        int nz_index; /* Fill matrix row of (i,j) */
+        row = center_index;
+        ptr = center_nz_index;
+
+        const Stencil& CenterStencil = getStencil(i_r);
+
+        offset = CenterStencil[StencilPosition::Center];
+        column = center_index;
+        value  = 0.25 * (h1 + h2) * (k1 + k2) * coeff_beta * std::fabs(detDF); /* beta_{i,j} */
+        update_CSR_COO_MatrixElement(matrix, ptr, offset, row, column, value);
+
+        offset = CenterStencil[StencilPosition::Left];
+        column = left_index;
+        value  = -coeff1 * arr; /* Left */
+        update_CSR_COO_MatrixElement(matrix, ptr, offset, row, column, value);
+
+        offset = CenterStencil[StencilPosition::Bottom];
+        column = bottom_index;
+        value  = -coeff3 * att; /* Bottom */
+        update_CSR_COO_MatrixElement(matrix, ptr, offset, row, column, value);
+
+        offset = CenterStencil[StencilPosition::Top];
+        column = top_index;
+        value  = -coeff4 * att; /* Top */
+        update_CSR_COO_MatrixElement(matrix, ptr, offset, row, column, value);
+
+        offset = CenterStencil[StencilPosition::Center];
+        column = center_index;
+        value  = (coeff1 + coeff2) * arr + (coeff3 + coeff4) * att; /* Center: (Left, Right, Bottom, Top) */
+        update_CSR_COO_MatrixElement(matrix, ptr, offset, row, column, value);
+
+        /* Fill matrix row of (i-1,j) */
+        /* From view the view of the across origin node, */ /* the directions are roatated by 180 degrees in the stencil! */
+        row = left_index;
+        ptr = left_nz_index;
+
+        const Stencil& LeftStencil = CenterStencil;
+
+        offset = LeftStencil[StencilPosition::Left];
+        column = center_index;
+        value  = -coeff1 * arr; /* Right -> Left*/
+        update_CSR_COO_MatrixElement(matrix, ptr, offset, row, column, value);
+
+        offset = LeftStencil[StencilPosition::Center];
+        column = left_index;
+        value  = +coeff1 * arr; /* Center: (Right) -> Center: (Left) */
+        update_CSR_COO_MatrixElement(matrix, ptr, offset, row, column, value);
+
+        /* Top Right -> Bottom Left: REMOVED DUE TO ARTIFICAL 7 POINT STENCIL */
+
+        /* Bottom Right -> Top Left: REMOVED DUE TO ARTIFICAL 7 POINT STENCIL */
+
+        /* Fill matrix row of (i,j-1) */
+        row = bottom_index;
+        ptr = bottom_nz_index;
+
+        const Stencil& BottomStencil = CenterStencil;
+
+        offset = BottomStencil[StencilPosition::Top];
+        column = center_index;
+        value  = -coeff3 * att; /* Top */
+        update_CSR_COO_MatrixElement(matrix, ptr, offset, row, column, value);
+
+        offset = BottomStencil[StencilPosition::Center];
+        column = bottom_index;
+        value  = +coeff3 * att; /* Center: (Top) */
+        update_CSR_COO_MatrixElement(matrix, ptr, offset, row, column, value);
+
+        /* TopLeft: REMOVED DUE TO ARTIFICAL 7 POINT STENCIL */
+
+        /* Fill matrix row of (i,j+1) */
+        row = top_index;
+        ptr = top_nz_index;
+
+        const Stencil& TopStencil = CenterStencil;
+
+        offset = TopStencil[StencilPosition::Bottom];
+        column = center_index;
+        value  = -coeff4 * att; /* Bottom */
+        update_CSR_COO_MatrixElement(matrix, ptr, offset, row, column, value);
+
+        offset = TopStencil[StencilPosition::Center];
+        column = top_index;
+        value  = +coeff4 * att; /* Center: (Bottom) */
+        update_CSR_COO_MatrixElement(matrix, ptr, offset, row, column, value);
+
+        /* BottomLeft: REMOVED DUE TO ARTIFICAL 7 POINT STENCIL */
+    }
+}
+
+template <concepts::DomainGeometry DomainGeometry>
+void SmootherGive<DomainGeometry>::nodeBuildInteriorBoundarySolverMatrix_i_r_1(int i_theta, const PolarGrid& grid,
+                                                                               bool DirBC_Interior,
+                                                                               InnerBoundaryMatrix& matrix, double arr,
+                                                                               double att, double art, double detDF,
+                                                                               double coeff_beta)
+{
+    using smoother_give::update_CSR_COO_MatrixElement;
+
+    assert(i_theta >= 0 && i_theta < grid.ntheta());
+
+    int ptr, offset;
+    int row, column;
+    double value;
+
+    const int i_r = 1;
+
+    const double h1 = grid.radialSpacing(i_r - 1);
+    const double h2 = grid.radialSpacing(i_r);
+    const double k1 = grid.angularSpacing(i_theta - 1);
+    const double k2 = grid.angularSpacing(i_theta);
+
+    const 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_index = i_theta;
+
+    /* Fill matrix row of (i-1,j) */
+    if (!DirBC_Interior) {
+        row = left_index;
+        ptr = getCircleAscIndex(i_r - 1, i_theta);
+
+        const Stencil& LeftStencil = getStencil(i_r - 1);
+
+        offset = LeftStencil[StencilPosition::Center];
+        column = left_index;
+        value  = coeff1 * arr; /* Center: (Right) */
+        update_CSR_COO_MatrixElement(matrix, ptr, offset, row, column, value);
+    }
+}
+
+template <concepts::DomainGeometry DomainGeometry>
+typename SmootherGive<DomainGeometry>::InnerBoundaryMatrix
+SmootherGive<DomainGeometry>::buildInteriorBoundarySolverMatrix()
+{
+    const PolarGrid& grid                         = Smoother<DomainGeometry>::grid_;
+    const LevelCache<DomainGeometry>& level_cache = Smoother<DomainGeometry>::level_cache_;
+    const bool DirBC_Interior                     = Smoother<DomainGeometry>::DirBC_Interior_;
+
+    const int i_r    = 0;
+    const int ntheta = grid.ntheta();
+
+#ifdef GMGPOLAR_USE_MUMPS
+    // Although the matrix is symmetric, we need to store all its entries, so we disable the symmetry.
+    const int nnz = getNonZeroCountCircleAsc(i_r);
+    SparseMatrixCOO<double> inner_boundary_solver_matrix(ntheta, ntheta, nnz);
+    inner_boundary_solver_matrix.is_symmetric(true);
+#else
+    // The stencils size for the inner boundary matrix is either 1 (Dirichlet BC) or 4 (across-origin discretization).
+    std::function<int(int)> nnz_per_row = [&](int i_theta) {
+        return DirBC_Interior ? 1 : 4;
+    };
+    SparseMatrixCSR<double> inner_boundary_solver_matrix(ntheta, ntheta, nnz_per_row);
+#endif
+
+    {
+        const int i_r  = 0;
+        const double r = grid.radius(i_r);
+        for (int i_theta = 0; i_theta < ntheta; i_theta++) {
+            {
+                const int global_index = grid.index(i_r, i_theta);
+                const double theta     = grid.theta(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);
+
+                nodeBuildInteriorBoundarySolverMatrix_i_r_0(i_theta, grid, DirBC_Interior, inner_boundary_solver_matrix,
+                                                            arr, att, art, detDF, coeff_beta);
+            }
+        }
+    }
+
+    {
+        const int i_r  = 1;
+        const double r = grid.radius(i_r);
+        for (int i_theta = 0; i_theta < ntheta; i_theta++) {
+            {
+                const int global_index = grid.index(i_r, i_theta);
+                const double theta     = grid.theta(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);
+
+                nodeBuildInteriorBoundarySolverMatrix_i_r_1(i_theta, grid, DirBC_Interior, inner_boundary_solver_matrix,
+                                                            arr, att, art, detDF, coeff_beta);
+            }
+        }
+    }
+
+    return inner_boundary_solver_matrix;
+}