Allow updating coefficients and reuse same solver instances. (Draft)

include/DirectSolver/DirectSolverGive/buildSolverMatrix.inl

@@ -794,29 +794,39 @@ typename DirectSolverGive<LevelCacheType>::SystemMatrix DirectSolverGive<LevelCa
 {
     using direct_solver_give::getNonZeroCountSolverMatrix;
     using direct_solver_give::getStencilSize;
-    using direct_solver_give::nodeBuildSolverMatrixGive;
-    using direct_solver_give::validateSolverMatrixIndexing;
 
-    const PolarGrid& grid             = DirectSolver<LevelCacheType>::grid_;
-    const LevelCacheType& level_cache = DirectSolver<LevelCacheType>::level_cache_;
-    const bool DirBC_Interior         = DirectSolver<LevelCacheType>::DirBC_Interior_;
-
-    assert(validateSolverMatrixIndexing(grid, DirBC_Interior) && "Solver matrix indexing is inconsistent");
-
-    const int n = grid.numberOfNodes();
+    const PolarGrid& grid     = DirectSolver<LevelCacheType>::grid_;
+    const bool DirBC_Interior = DirectSolver<LevelCacheType>::DirBC_Interior_;
+    const int n               = grid.numberOfNodes();
 
 #ifdef GMGPOLAR_USE_MUMPS
     const int nnz = getNonZeroCountSolverMatrix(grid, DirBC_Interior);
-    SparseMatrixCOO<double> solver_matrix(n, n, nnz);
+    SystemMatrix solver_matrix(n, n, nnz);
     solver_matrix.is_symmetric(true);
 #else
     std::function<int(int)> nnz_per_row = [&](int global_index) {
         return getStencilSize(global_index, grid, DirBC_Interior);
     };
-
-    SparseMatrixCSR<double> solver_matrix(n, n, nnz_per_row);
+    SystemMatrix solver_matrix(n, n, nnz_per_row);
 #endif
 
+    fillSolverMatrix(solver_matrix);
+
+    return solver_matrix;
+}
+
+template <class LevelCacheType>
+void DirectSolverGive<LevelCacheType>::fillSolverMatrix(SystemMatrix& solver_matrix)
+{
+    using direct_solver_give::nodeBuildSolverMatrixGive;
+    using direct_solver_give::validateSolverMatrixIndexing;
+
+    const PolarGrid& grid             = DirectSolver<LevelCacheType>::grid_;
+    const LevelCacheType& level_cache = DirectSolver<LevelCacheType>::level_cache_;
+    const bool DirBC_Interior         = DirectSolver<LevelCacheType>::DirBC_Interior_;
+
+    assert(validateSolverMatrixIndexing(grid, DirBC_Interior) && "Solver matrix indexing is inconsistent");
+
     /* ---------------- */
     /* Circular section */
     /* ---------------- */
@@ -872,6 +882,4 @@ typename DirectSolverGive<LevelCacheType>::SystemMatrix DirectSolverGive<LevelCa
             });
         Kokkos::fence();
     }
-
-    return solver_matrix;
 }

include/DirectSolver/DirectSolverGive/directSolverGive.h

@@ -13,6 +13,7 @@ class DirectSolverGive : public DirectSolver<LevelCacheType>
 
     // Note: The rhs (right-hand side) vector gets overwritten during the solution process.
     void solveInPlace(Vector<double> solution) override;
+    void updateValues() override;
 
 private:
 #ifdef GMGPOLAR_USE_MUMPS
@@ -21,11 +22,8 @@ class DirectSolverGive : public DirectSolver<LevelCacheType>
 #else
     using SystemMatrix = SparseMatrixCSR<double>;
     using SystemSolver = SparseLUSolver<double>;
-    // Stored only for the in-house solver (CSR).
-    SystemMatrix system_matrix_;
 #endif
-
-    // Solver object (owns matrix if MUMPS, references if in-house solver).
+    SystemMatrix system_matrix_;
     SystemSolver system_solver_;
 
 public:
@@ -42,6 +40,13 @@ class DirectSolverGive : public DirectSolver<LevelCacheType>
     void applySymmetryShift(Vector<double> rhs) const;
     void applySymmetryShiftInnerBoundary(Vector<double> x) const;
     void applySymmetryShiftOuterBoundary(Vector<double> x) const;
+
+private:
+    // Shared kernel: runs the colored circular/radial passes to fill
+    // solver_matrix's entries from the current grid + LevelCache data.
+    // Used by buildSolverMatrix() (fresh storage) and updateValues()
+    // (existing storage, values only).
+    void fillSolverMatrix(SystemMatrix& solver_matrix);
 };
 
 #include "applySymmetryShift.inl"

include/DirectSolver/DirectSolverGive/directSolverGive.inl

@@ -4,12 +4,8 @@ template <class LevelCacheType>
 DirectSolverGive<LevelCacheType>::DirectSolverGive(const PolarGrid& grid, const LevelCacheType& level_cache,
                                                    bool DirBC_Interior)
     : DirectSolver<LevelCacheType>(grid, level_cache, DirBC_Interior)
-#ifdef GMGPOLAR_USE_MUMPS
-    , system_solver_(buildSolverMatrix())
-#else
     , system_matrix_(buildSolverMatrix())
     , system_solver_(system_matrix_)
-#endif
 {
 }
 
@@ -25,3 +21,16 @@ void DirectSolverGive<LevelCacheType>::solveInPlace(Vector<double> solution)
     // Solves the adjusted system symmetric(matrixA) * solution = rhs using the MUMPS solver.
     system_solver_.solveInPlace(solution);
 }
+
+template <class LevelCacheType>
+void DirectSolverGive<LevelCacheType>::updateValues()
+{
+    // Overwrite system_matrix_'s numeric entries in place. Sparsity pattern
+    // is unchanged (same grid, same stencil coloring), so no reallocation
+    // occurs here.
+    fillSolverMatrix(system_matrix_);
+
+    // Refactorize, reusing whatever one-time analysis/ordering work each
+    // backend already did (RCM ordering / MUMPS analysis phase).
+    system_solver_.updateValues(system_matrix_);
+}

include/DirectSolver/DirectSolverTake/buildSolverMatrix.inl

@@ -475,29 +475,40 @@ typename DirectSolverTake<LevelCacheType>::SystemMatrix DirectSolverTake<LevelCa
 {
     using direct_solver_take::getNonZeroCountSolverMatrix;
     using direct_solver_take::getStencilSize;
-    using direct_solver_take::nodeBuildSolverMatrixTake;
-    using direct_solver_take::validateSolverMatrixIndexing;
 
-    const PolarGrid& grid             = DirectSolver<LevelCacheType>::grid_;
-    const LevelCacheType& level_cache = DirectSolver<LevelCacheType>::level_cache_;
-    const bool DirBC_Interior         = DirectSolver<LevelCacheType>::DirBC_Interior_;
-
-    assert(validateSolverMatrixIndexing(grid, DirBC_Interior) && "Solver matrix indexing is inconsistent");
+    const PolarGrid& grid     = DirectSolver<LevelCacheType>::grid_;
+    const bool DirBC_Interior = DirectSolver<LevelCacheType>::DirBC_Interior_;
 
     const int n = grid.numberOfNodes();
 
 #ifdef GMGPOLAR_USE_MUMPS
     const int nnz = getNonZeroCountSolverMatrix(grid, DirBC_Interior);
-    SparseMatrixCOO<double> solver_matrix(n, n, nnz);
+    SystemMatrix solver_matrix(n, n, nnz);
     solver_matrix.is_symmetric(true);
 #else
     std::function<int(int)> nnz_per_row = [&](int global_index) {
         return getStencilSize(global_index, grid, DirBC_Interior);
     };
-
-    SparseMatrixCSR<double> solver_matrix(n, n, nnz_per_row);
+    SystemMatrix solver_matrix(n, n, nnz_per_row);
 #endif
 
+    fillSolverMatrix(solver_matrix); // allocation done above, this only fills entries
+
+    return solver_matrix;
+}
+
+template <class LevelCacheType>
+void DirectSolverTake<LevelCacheType>::fillSolverMatrix(SystemMatrix& solver_matrix)
+{
+    using direct_solver_take::nodeBuildSolverMatrixTake;
+    using direct_solver_take::validateSolverMatrixIndexing;
+
+    const PolarGrid& grid             = DirectSolver<LevelCacheType>::grid_;
+    const LevelCacheType& level_cache = DirectSolver<LevelCacheType>::level_cache_;
+    const bool DirBC_Interior         = DirectSolver<LevelCacheType>::DirBC_Interior_;
+
+    assert(validateSolverMatrixIndexing(grid, DirBC_Interior) && "Solver matrix indexing is inconsistent");
+
     assert(level_cache.cacheDensityProfileCoefficients());
     assert(level_cache.cacheDomainGeometry());
 
@@ -537,6 +548,4 @@ typename DirectSolverTake<LevelCacheType>::SystemMatrix DirectSolverTake<LevelCa
         });
 
     Kokkos::fence();
-
-    return solver_matrix;
-}
+}

include/DirectSolver/DirectSolverTake/directSolverTake.h

@@ -13,6 +13,7 @@ class DirectSolverTake : public DirectSolver<LevelCacheType>
 
     // Note: The rhs (right-hand side) vector gets overwritten during the solution process.
     void solveInPlace(Vector<double> solution) override;
+    void updateValues() override;
 
 private:
 #ifdef GMGPOLAR_USE_MUMPS
@@ -21,11 +22,8 @@ class DirectSolverTake : public DirectSolver<LevelCacheType>
 #else
     using SystemMatrix = SparseMatrixCSR<double>;
     using SystemSolver = SparseLUSolver<double>;
-    // Stored only for the in-house solver (CSR).
-    SystemMatrix system_matrix_;
 #endif
-
-    // Solver object (owns matrix if MUMPS, references if in-house solver).
+    SystemMatrix system_matrix_;
     SystemSolver system_solver_;
 
 public:
@@ -42,6 +40,12 @@ class DirectSolverTake : public DirectSolver<LevelCacheType>
     void applySymmetryShift(Vector<double> rhs) const;
     void applySymmetryShiftInnerBoundary(Vector<double> rhs) const;
     void applySymmetryShiftOuterBoundary(Vector<double> rhs) const;
+
+private:
+    // Shared kernel: fills solver_matrix's entries from the current grid +
+    // LevelCache data. Used by buildSolverMatrix() (on freshly allocated
+    // storage) and updateValues() (on existing storage, values only).
+    void fillSolverMatrix(SystemMatrix& solver_matrix);
 };
 
 #include "applySymmetryShift.inl"

include/DirectSolver/DirectSolverTake/directSolverTake.inl

@@ -1,27 +1,55 @@
 #pragma once
 
+// template <class LevelCacheType>
+// DirectSolverTake<LevelCacheType>::DirectSolverTake(const PolarGrid& grid, const LevelCacheType& level_cache,
+//                                                    bool DirBC_Interior)
+//     : DirectSolver<LevelCacheType>(grid, level_cache, DirBC_Interior)
+// #ifdef GMGPOLAR_USE_MUMPS
+//     , system_solver_(buildSolverMatrix())
+// #else
+//     , system_matrix_(buildSolverMatrix())
+//     , system_solver_(system_matrix_)
+// #endif
+// {
+// }
+
+// template <class LevelCacheType>
+// void DirectSolverTake<LevelCacheType>::solveInPlace(Vector<double> solution)
+// {
+//     // Adjusts the right-hand side vector to account for symmetry corrections.
+//     // This transforms the system matrixA * solution = rhs into the equivalent system:
+//     // symmetric_DBc(matrixA) * solution = rhs - applySymmetryShift(rhs).
+//     // The correction modifies the rhs to account for the influence of the Dirichlet boundary conditions,
+//     // ensuring that the solution at the boundary is correctly adjusted and maintains the required symmetry.
+//     applySymmetryShift(solution);
+//     // Solves the adjusted system symmetric(matrixA) * solution = rhs using the MUMPS solver.
+//     system_solver_.solveInPlace(solution);
+// }
+
 template <class LevelCacheType>
 DirectSolverTake<LevelCacheType>::DirectSolverTake(const PolarGrid& grid, const LevelCacheType& level_cache,
                                                    bool DirBC_Interior)
     : DirectSolver<LevelCacheType>(grid, level_cache, DirBC_Interior)
-#ifdef GMGPOLAR_USE_MUMPS
-    , system_solver_(buildSolverMatrix())
-#else
     , system_matrix_(buildSolverMatrix())
     , system_solver_(system_matrix_)
-#endif
 {
 }
 
 template <class LevelCacheType>
 void DirectSolverTake<LevelCacheType>::solveInPlace(Vector<double> solution)
 {
-    // Adjusts the right-hand side vector to account for symmetry corrections.
-    // This transforms the system matrixA * solution = rhs into the equivalent system:
-    // symmetric_DBc(matrixA) * solution = rhs - applySymmetryShift(rhs).
-    // The correction modifies the rhs to account for the influence of the Dirichlet boundary conditions,
-    // ensuring that the solution at the boundary is correctly adjusted and maintains the required symmetry.
     applySymmetryShift(solution);
-    // Solves the adjusted system symmetric(matrixA) * solution = rhs using the MUMPS solver.
     system_solver_.solveInPlace(solution);
 }
+
+template <class LevelCacheType>
+void DirectSolverTake<LevelCacheType>::updateValues()
+{
+    // Refill system_matrix_'s numeric entries in place; sparsity pattern
+    // (grid + DirBC_Interior) is unchanged, so no reallocation happens here.
+    fillSolverMatrix(system_matrix_);
+
+    // Refactorize, reusing whatever one-time analysis/ordering work each
+    // backend already did (RCM ordering / MUMPS analysis phase).
+    system_solver_.updateValues(system_matrix_);
+}

include/DirectSolver/directSolver.h

@@ -38,6 +38,12 @@ class DirectSolver
     // Note: The rhs (right-hand side) vector gets overwritten during the solution process.
     virtual void solveInPlace(Vector<double> solution) = 0;
 
+    // Recomputes the solver matrix entries from the current LevelCache values
+    // (arr, att, art, detDF, coeff_beta) and refreshes the factorization.
+    // No reallocation: sparsity pattern and matrix storage are unchanged,
+    // only the numeric values are overwritten.
+    virtual void updateValues() = 0;
+
 protected:
     const PolarGrid& grid_;
     const LevelCacheType& level_cache_;

include/GMGPolar/gmgpolar.h

@@ -33,8 +33,8 @@ class GMGPolar : public IGMGPolar
     GMGPolar(const PolarGrid& grid, const DomainGeometry& domain_geometry,
              const DensityProfileCoefficients& density_profile_coefficients)
         : IGMGPolar(grid)
-        , domain_geometry_(domain_geometry)
-        , density_profile_coefficients_(density_profile_coefficients)
+        , domain_geometry_(&domain_geometry)
+        , density_profile_coefficients_(&density_profile_coefficients)
         , exact_solution_(nullptr)
         // Level management and internal solver data
         , number_of_levels_(0)
@@ -54,6 +54,9 @@ class GMGPolar : public IGMGPolar
     // Finalize solver setup (allocate data, build operators, etc.).
     void setup();
 
+    void updateOperatorValues(const DomainGeometry* domain_geometry                          = nullptr,
+                              const DensityProfileCoefficients* density_profile_coefficients = nullptr);
+
     // Solve system with given boundary conditions and source term.
     // Multiple solves with different inputs are supported.
     template <concepts::BoundaryConditions BoundaryConditions, concepts::SourceTerm SourceTerm>
@@ -89,8 +92,8 @@ class GMGPolar : public IGMGPolar
     /* ------------------------------------ */
     /* Grid Configuration & Input Functions */
     /* ------------------------------------ */
-    const DomainGeometry& domain_geometry_;
-    const DensityProfileCoefficients& density_profile_coefficients_;
+    const DomainGeometry* domain_geometry_;
+    const DensityProfileCoefficients* density_profile_coefficients_;
     const ExactSolution* exact_solution_; // Optional exact solution for validation
 
     /* ---------------- */
@@ -151,6 +154,8 @@ class GMGPolar : public IGMGPolar
                      const SourceTerm& source_term);
     void discretize_rhs_f(const LevelType& level, Vector<double> rhs_f);
 
+    void recomputeOperatorValues();
+
     /* --------------- */
     /* Solve Functions */
     void applyExtrapolation(int current_level, Vector<double> fine_values, ConstVector<double> coarse_values);

include/GMGPolar/setup.h

@@ -41,14 +41,14 @@ void GMGPolar<DomainGeometry, DensityProfileCoefficients>::setup()
     levels_.reserve(number_of_levels_);
 
     auto finest_levelCache = std::make_unique<LevelCache<DomainGeometry, DensityProfileCoefficients>>(
-        *finest_grid, density_profile_coefficients_, domain_geometry_, cache_density_profile_coefficients_,
+        *finest_grid, *density_profile_coefficients_, *domain_geometry_, cache_density_profile_coefficients_,
         cache_domain_geometry_, 0);
     levels_.emplace_back(0, std::move(finest_grid), std::move(finest_levelCache), extrapolation_, FMG_, PCG_FMG_);
 
     for (int level_depth = 1; level_depth < number_of_levels_; level_depth++) {
         auto current_grid       = std::make_unique<PolarGrid>(coarseningGrid(levels_[level_depth - 1].grid()));
         auto current_levelCache = std::make_unique<LevelCache<DomainGeometry, DensityProfileCoefficients>>(
-            *current_grid, density_profile_coefficients_, domain_geometry_, cache_density_profile_coefficients_,
+            *current_grid, *density_profile_coefficients_, *domain_geometry_, cache_density_profile_coefficients_,
             cache_domain_geometry_, level_depth);
         levels_.emplace_back(level_depth, std::move(current_grid), std::move(current_levelCache), extrapolation_, FMG_,
                              PCG_FMG_);
@@ -253,7 +253,7 @@ void GMGPolar<DomainGeometry, DensityProfileCoefficients>::discretize_rhs_f(cons
     else {
         /* DomainGeometry is not cached */
         // Local copy is required to avoid copying the class
-        const DomainGeometry& domain_geometry = domain_geometry_;
+        const DomainGeometry& domain_geometry = *domain_geometry_;
 
         // ---------------------------------------------- //
         // Discretize rhs values (circular index section) //
@@ -385,6 +385,29 @@ void GMGPolar<DomainGeometry, DensityProfileCoefficients>::build_rhs_f(const Lev
     Kokkos::fence();
 }
 
+template <concepts::DomainGeometry DomainGeometry, concepts::DensityProfileCoefficients DensityProfileCoefficients>
+void GMGPolar<DomainGeometry, DensityProfileCoefficients>::updateOperatorValues(
+    const DomainGeometry* domain_geometry, const DensityProfileCoefficients* density_profile_coefficients)
+{
+    if (domain_geometry) {
+        domain_geometry_ = domain_geometry;
+    }
+    if (density_profile_coefficients) {
+        density_profile_coefficients_ = density_profile_coefficients;
+    }
+
+    recomputeOperatorValues();
+}
+
+template <concepts::DomainGeometry DomainGeometry, concepts::DensityProfileCoefficients DensityProfileCoefficients>
+void GMGPolar<DomainGeometry, DensityProfileCoefficients>::recomputeOperatorValues()
+
+{
+    for (int level_depth = 0; level_depth < number_of_levels_; ++level_depth) {
+        levels_[level_depth].updateOperatorValues(*domain_geometry_, *density_profile_coefficients_);
+    }
+}
+
 template <concepts::DomainGeometry DomainGeometry, concepts::DensityProfileCoefficients DensityProfileCoefficients>
 void GMGPolar<DomainGeometry, DensityProfileCoefficients>::printSettings(const PolarGrid& finest_grid,
                                                                          const PolarGrid& coarsest_grid) const