From 308583544925dd9b96c7e1901a6d3f846ad5b430 Mon Sep 17 00:00:00 2001 From: julianlitz Date: Tue, 11 Aug 2026 22:54:20 +0200 Subject: [PATCH 1/2] First Draft --- .../DirectSolverGive/buildSolverMatrix.inl | 36 +++++++----- .../DirectSolverGive/directSolverGive.h | 13 +++-- .../DirectSolverGive/directSolverGive.inl | 17 ++++-- .../DirectSolverTake/buildSolverMatrix.inl | 35 +++++++----- .../DirectSolverTake/directSolverTake.h | 7 +++ .../DirectSolverTake/directSolverTake.inl | 48 ++++++++++++---- include/DirectSolver/directSolver.h | 6 ++ include/GMGPolar/gmgpolar.h | 13 +++-- include/GMGPolar/setup.h | 29 +++++++++- include/GMGPolar/utils.h | 6 +- include/Level/level.h | 26 ++++++--- include/Level/level.inl | 40 +++++++++++-- include/Level/levelCache.inl | 26 +++++++++ .../LinearAlgebra/Solvers/coo_mumps_solver.h | 56 +++++++++++++++++++ include/LinearAlgebra/Solvers/csr_lu_solver.h | 30 ++++++++++ scripts/compile.sh | 2 +- 16 files changed, 324 insertions(+), 66 deletions(-) diff --git a/include/DirectSolver/DirectSolverGive/buildSolverMatrix.inl b/include/DirectSolver/DirectSolverGive/buildSolverMatrix.inl index a1c02bea0..83817448e 100644 --- a/include/DirectSolver/DirectSolverGive/buildSolverMatrix.inl +++ b/include/DirectSolver/DirectSolverGive/buildSolverMatrix.inl @@ -794,29 +794,39 @@ typename DirectSolverGive::SystemMatrix DirectSolverGive::grid_; - const LevelCacheType& level_cache = DirectSolver::level_cache_; - const bool DirBC_Interior = DirectSolver::DirBC_Interior_; - - assert(validateSolverMatrixIndexing(grid, DirBC_Interior) && "Solver matrix indexing is inconsistent"); - - const int n = grid.numberOfNodes(); + const PolarGrid& grid = DirectSolver::grid_; + const bool DirBC_Interior = DirectSolver::DirBC_Interior_; + const int n = grid.numberOfNodes(); #ifdef GMGPOLAR_USE_MUMPS const int nnz = getNonZeroCountSolverMatrix(grid, DirBC_Interior); - SparseMatrixCOO solver_matrix(n, n, nnz); + SystemMatrix solver_matrix(n, n, nnz); solver_matrix.is_symmetric(true); #else std::function nnz_per_row = [&](int global_index) { return getStencilSize(global_index, grid, DirBC_Interior); }; - - SparseMatrixCSR solver_matrix(n, n, nnz_per_row); + SystemMatrix solver_matrix(n, n, nnz_per_row); #endif + fillSolverMatrix(solver_matrix); + + return solver_matrix; +} + +template +void DirectSolverGive::fillSolverMatrix(SystemMatrix& solver_matrix) +{ + using direct_solver_give::nodeBuildSolverMatrixGive; + using direct_solver_give::validateSolverMatrixIndexing; + + const PolarGrid& grid = DirectSolver::grid_; + const LevelCacheType& level_cache = DirectSolver::level_cache_; + const bool DirBC_Interior = DirectSolver::DirBC_Interior_; + + assert(validateSolverMatrixIndexing(grid, DirBC_Interior) && "Solver matrix indexing is inconsistent"); + /* ---------------- */ /* Circular section */ /* ---------------- */ @@ -872,6 +882,4 @@ typename DirectSolverGive::SystemMatrix DirectSolverGive // Note: The rhs (right-hand side) vector gets overwritten during the solution process. void solveInPlace(Vector solution) override; + void updateValues() override; private: #ifdef GMGPOLAR_USE_MUMPS @@ -21,11 +22,8 @@ class DirectSolverGive : public DirectSolver #else using SystemMatrix = SparseMatrixCSR; using SystemSolver = SparseLUSolver; - // 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 void applySymmetryShift(Vector rhs) const; void applySymmetryShiftInnerBoundary(Vector x) const; void applySymmetryShiftOuterBoundary(Vector 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" diff --git a/include/DirectSolver/DirectSolverGive/directSolverGive.inl b/include/DirectSolver/DirectSolverGive/directSolverGive.inl index 6d0a74926..44d7eb7bb 100644 --- a/include/DirectSolver/DirectSolverGive/directSolverGive.inl +++ b/include/DirectSolver/DirectSolverGive/directSolverGive.inl @@ -4,12 +4,8 @@ template DirectSolverGive::DirectSolverGive(const PolarGrid& grid, const LevelCacheType& level_cache, bool DirBC_Interior) : DirectSolver(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::solveInPlace(Vector solution) // Solves the adjusted system symmetric(matrixA) * solution = rhs using the MUMPS solver. system_solver_.solveInPlace(solution); } + +template +void DirectSolverGive::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_); +} \ No newline at end of file diff --git a/include/DirectSolver/DirectSolverTake/buildSolverMatrix.inl b/include/DirectSolver/DirectSolverTake/buildSolverMatrix.inl index 3f04c6104..ae6fb4446 100644 --- a/include/DirectSolver/DirectSolverTake/buildSolverMatrix.inl +++ b/include/DirectSolver/DirectSolverTake/buildSolverMatrix.inl @@ -475,29 +475,40 @@ typename DirectSolverTake::SystemMatrix DirectSolverTake::grid_; - const LevelCacheType& level_cache = DirectSolver::level_cache_; - const bool DirBC_Interior = DirectSolver::DirBC_Interior_; - - assert(validateSolverMatrixIndexing(grid, DirBC_Interior) && "Solver matrix indexing is inconsistent"); + const PolarGrid& grid = DirectSolver::grid_; + const bool DirBC_Interior = DirectSolver::DirBC_Interior_; const int n = grid.numberOfNodes(); #ifdef GMGPOLAR_USE_MUMPS const int nnz = getNonZeroCountSolverMatrix(grid, DirBC_Interior); - SparseMatrixCOO solver_matrix(n, n, nnz); + SystemMatrix solver_matrix(n, n, nnz); solver_matrix.is_symmetric(true); #else std::function nnz_per_row = [&](int global_index) { return getStencilSize(global_index, grid, DirBC_Interior); }; - - SparseMatrixCSR 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 +void DirectSolverTake::fillSolverMatrix(SystemMatrix& solver_matrix) +{ + using direct_solver_take::nodeBuildSolverMatrixTake; + using direct_solver_take::validateSolverMatrixIndexing; + + const PolarGrid& grid = DirectSolver::grid_; + const LevelCacheType& level_cache = DirectSolver::level_cache_; + const bool DirBC_Interior = DirectSolver::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::SystemMatrix DirectSolverTake // Note: The rhs (right-hand side) vector gets overwritten during the solution process. void solveInPlace(Vector solution) override; + void updateValues() override; private: #ifdef GMGPOLAR_USE_MUMPS @@ -42,6 +43,12 @@ class DirectSolverTake : public DirectSolver void applySymmetryShift(Vector rhs) const; void applySymmetryShiftInnerBoundary(Vector rhs) const; void applySymmetryShiftOuterBoundary(Vector 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" diff --git a/include/DirectSolver/DirectSolverTake/directSolverTake.inl b/include/DirectSolver/DirectSolverTake/directSolverTake.inl index f7efac3c8..99bbcb724 100644 --- a/include/DirectSolver/DirectSolverTake/directSolverTake.inl +++ b/include/DirectSolver/DirectSolverTake/directSolverTake.inl @@ -1,27 +1,55 @@ #pragma once +// template +// DirectSolverTake::DirectSolverTake(const PolarGrid& grid, const LevelCacheType& level_cache, +// bool DirBC_Interior) +// : DirectSolver(grid, level_cache, DirBC_Interior) +// #ifdef GMGPOLAR_USE_MUMPS +// , system_solver_(buildSolverMatrix()) +// #else +// , system_matrix_(buildSolverMatrix()) +// , system_solver_(system_matrix_) +// #endif +// { +// } + +// template +// void DirectSolverTake::solveInPlace(Vector 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 DirectSolverTake::DirectSolverTake(const PolarGrid& grid, const LevelCacheType& level_cache, bool DirBC_Interior) : DirectSolver(grid, level_cache, DirBC_Interior) -#ifdef GMGPOLAR_USE_MUMPS - , system_solver_(buildSolverMatrix()) -#else , system_matrix_(buildSolverMatrix()) , system_solver_(system_matrix_) -#endif { } template void DirectSolverTake::solveInPlace(Vector 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 +void DirectSolverTake::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_); +} \ No newline at end of file diff --git a/include/DirectSolver/directSolver.h b/include/DirectSolver/directSolver.h index ce43fa129..1fca471ba 100644 --- a/include/DirectSolver/directSolver.h +++ b/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 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_; diff --git a/include/GMGPolar/gmgpolar.h b/include/GMGPolar/gmgpolar.h index daed475ac..e6b7bc279 100644 --- a/include/GMGPolar/gmgpolar.h +++ b/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 @@ -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 rhs_f); + void recomputeOperatorValues(); + /* --------------- */ /* Solve Functions */ void applyExtrapolation(int current_level, Vector fine_values, ConstVector coarse_values); diff --git a/include/GMGPolar/setup.h b/include/GMGPolar/setup.h index 250ef3979..d359d8653 100644 --- a/include/GMGPolar/setup.h +++ b/include/GMGPolar/setup.h @@ -41,14 +41,14 @@ void GMGPolar::setup() levels_.reserve(number_of_levels_); auto finest_levelCache = std::make_unique>( - *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(coarseningGrid(levels_[level_depth - 1].grid())); auto current_levelCache = std::make_unique>( - *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::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::build_rhs_f(const Lev Kokkos::fence(); } +template +void GMGPolar::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 +void GMGPolar::recomputeOperatorValues() + +{ + for (int level_depth = 0; level_depth < number_of_levels_; ++level_depth) { + levels_[level_depth].updateOperatorValues(*domain_geometry_, *density_profile_coefficients_); + } +} + template void GMGPolar::printSettings(const PolarGrid& finest_grid, const PolarGrid& coarsest_grid) const diff --git a/include/GMGPolar/utils.h b/include/GMGPolar/utils.h index 79f147098..491fd4b18 100644 --- a/include/GMGPolar/utils.h +++ b/include/GMGPolar/utils.h @@ -209,7 +209,7 @@ void GMGPolar::writeToVTK(const std: << "\n"; Vector Fx("Fx", grid.numberOfNodes()); Vector Fy("Fy", grid.numberOfNodes()); - const DomainGeometry& domain_geometry = domain_geometry_; + const DomainGeometry& domain_geometry = *domain_geometry_; Kokkos::parallel_for( "collect Fx,Fy", Kokkos::RangePolicy(0, grid.numberOfNodes()), KOKKOS_LAMBDA(int index) { @@ -280,13 +280,15 @@ void GMGPolar::writeToVTK(const std: // Write points file << "\n" << "\n"; + + const DomainGeometry& domain_geometry = *domain_geometry_; int i_r, i_theta; double r, theta; for (int index = 0; index < grid.numberOfNodes(); index++) { grid.multiIndex(index, i_r, i_theta); r = grid.radius(i_r); theta = grid.theta(i_theta); - file << domain_geometry_.Fx(r, theta) << " " << domain_geometry_.Fy(r, theta) << " " << 0 << "\n"; + file << domain_geometry.Fx(r, theta) << " " << domain_geometry.Fy(r, theta) << " " << 0 << "\n"; } file << "\n" diff --git a/include/Level/level.h b/include/Level/level.h index 101a37f98..39186d533 100644 --- a/include/Level/level.h +++ b/include/Level/level.h @@ -73,8 +73,14 @@ class Level // ----------- // // Constructor // explicit Level(const int level_depth, std::unique_ptr grid, - std::unique_ptr level_cache, const ExtrapolationType extrapolation, - const bool FMG, const bool PCG_FMG = false); + std::unique_ptr level_cache, const ExtrapolationType extrapolation, const bool FMG, + const bool PCG_FMG = false); + + // Pushes new geometry/coefficients into the LevelCache (overwriting cached + // values in place) and refreshes any multigrid operators that cache + // derived data. No reallocation. + void updateOperatorValues(const DomainGeometry& domain_geometry, + const DensityProfileCoefficients& density_profile_coefficients); // ---------------- // // Getter Functions // @@ -117,7 +123,10 @@ class Level private: const int level_depth_; std::unique_ptr grid_; - std::unique_ptr level_cache_; + std::unique_ptr level_cache_; + + // Refreshes cached operator data (matrix entries, factorizations, etc.) derived from the current LevelCache values + void updateMultigridOperators(); std::unique_ptr> op_directSolver_; std::unique_ptr> op_residual_; @@ -138,6 +147,9 @@ class LevelCache const DomainGeometry& domain_geometry, const bool cache_density_profile_coefficients, const bool cache_domain_geometry, const int level_depth); + void updateValues(const PolarGrid& grid, const DomainGeometry& domain_geometry, + const DensityProfileCoefficients& density_profile_coefficients); + const DomainGeometry& domainGeometry() const; const DensityProfileCoefficients& densityProfileCoefficients() const; @@ -157,8 +169,8 @@ class LevelCache { const int i_r_glob = i_r << level_depth_; const int i_theta_glob = i_theta << level_depth_; - coeff_beta = cache_density_profile_coefficients_ ? coeff_beta_[global_index] - : density_profile_coefficients_.beta(i_r_glob, i_theta_glob); + coeff_beta = cache_density_profile_coefficients_ ? coeff_beta_[global_index] + : density_profile_coefficients_.beta(i_r_glob, i_theta_glob); if (cache_domain_geometry_) { arr = arr_[global_index]; @@ -177,8 +189,8 @@ class LevelCache private: const int level_depth_; - const DomainGeometry domain_geometry_; - const DensityProfileCoefficients density_profile_coefficients_; + DomainGeometry domain_geometry_; + DensityProfileCoefficients density_profile_coefficients_; bool cache_density_profile_coefficients_; // cache alpha(r, theta), beta(r, theta) Vector coeff_alpha_; diff --git a/include/Level/level.inl b/include/Level/level.inl index 451487609..4176e22ba 100644 --- a/include/Level/level.inl +++ b/include/Level/level.inl @@ -3,10 +3,10 @@ // ----------- // // Constructor // template -Level::Level( - const int level_depth, std::unique_ptr grid, - std::unique_ptr> level_cache, - const ExtrapolationType extrapolation, const bool FMG, const bool PCG_FMG) +Level::Level(const int level_depth, std::unique_ptr grid, + std::unique_ptr level_cache, + const ExtrapolationType extrapolation, const bool FMG, + const bool PCG_FMG) : level_depth_(level_depth) , grid_(std::move(grid)) , level_cache_(std::move(level_cache)) @@ -19,6 +19,38 @@ Level::Level( { } +template +void Level::updateOperatorValues( + const DomainGeometry& domain_geometry, const DensityProfileCoefficients& density_profile_coefficients) +{ + // Overwrite cached geometry/coefficient values in place. + level_cache_->updateValues(*grid_, domain_geometry, density_profile_coefficients); + + // Refresh any multigrid operators built on top of those cached values. + updateMultigridOperators(); +} + +template +void Level::updateMultigridOperators() +{ + // TODO: Once DirectSolver / Residual / Smoother / ExtrapolatedSmoother expose + // an updateValues()-style method, call it here conditionally, mirroring the + // initialize*() pattern already used elsewhere in this class, e.g.: + // + // if (op_residual_) { + // op_residual_->updateValues(); + // } + if (op_directSolver_) { + op_directSolver_->updateValues(); + } + // if (op_smoother_) { + // op_smoother_->updateValues(); + // } + // if (op_extrapolated_smoother_) { + // op_extrapolated_smoother_->updateValues(); + // } +} + // ---------------- // // Getter Functions // template diff --git a/include/Level/levelCache.inl b/include/Level/levelCache.inl index 20a66e1d3..035b5d68a 100644 --- a/include/Level/levelCache.inl +++ b/include/Level/levelCache.inl @@ -2,6 +2,7 @@ namespace level_cache_helpers { + template static void cache_density_profile_coefficients(const PolarGrid& grid, const DensityProfileCoefficients& density_profile_coefficients, @@ -120,6 +121,31 @@ LevelCache::LevelCache( Kokkos::fence(); } +template +void LevelCache::updateValues( + const PolarGrid& grid, const DomainGeometry& domain_geometry, + const DensityProfileCoefficients& density_profile_coefficients) +{ + // Store the new geometry/coefficients (cheap object copy, not the expensive part). + domain_geometry_ = domain_geometry; + density_profile_coefficients_ = density_profile_coefficients; + + // Re-run the same caching kernels used at construction time. coeff_alpha_, + // coeff_beta_, arr_, att_, art_, detDF_ are already sized correctly and are + // simply overwritten in place - no Kokkos::View reallocation occurs here. + if (cache_density_profile_coefficients_) { + level_cache_helpers::cache_density_profile_coefficients(grid, density_profile_coefficients_, coeff_alpha_, + coeff_beta_, cache_domain_geometry_, level_depth_); + } + + if (cache_domain_geometry_) { + level_cache_helpers::cache_domain_geometry(grid, density_profile_coefficients_, domain_geometry_, arr_, att_, + art_, detDF_, level_depth_); + } + + Kokkos::fence(); +} + template const DensityProfileCoefficients& LevelCache::densityProfileCoefficients() const diff --git a/include/LinearAlgebra/Solvers/coo_mumps_solver.h b/include/LinearAlgebra/Solvers/coo_mumps_solver.h index c748aed39..e218b73a8 100644 --- a/include/LinearAlgebra/Solvers/coo_mumps_solver.h +++ b/include/LinearAlgebra/Solvers/coo_mumps_solver.h @@ -55,6 +55,62 @@ class CooMumpsSolver // rhs is overwritten in-place with the solution on return. void solveInPlace(Vector& rhs); + /** + * @brief Refactorize using a matrix with the same sparsity pattern as the + * one originally used to construct this solver. + * + * Reuses the existing MUMPS analysis/ordering phase (job 1) and only + * reruns the numeric factorization (job 2). The matrix must have the same + * non-zero structure - same number of entries in the same order - as the + * matrix originally passed to the constructor. + * + * Takes matrix by const reference (unlike the constructor, which takes + * ownership) so the caller can keep refilling and reusing their own + * matrix storage across repeated updates. + */ + template + void updateValues(const SparseMatrixCOO& matrix) + { + auto matrix_host = matrix.template mirror_view_and_copy(); + + if (matrix_host.is_symmetric()) { + matrix_ = extractUpperTriangle(matrix_host); + } + else { + matrix_ = std::move(matrix_host); + } + + assert(matrix_.non_zero_size() == mumps_solver_.nz && + "Matrix structure must match the originally factorized matrix"); + + // MUMPS uses 1-based indexing. + for (int i = 0; i < matrix_.non_zero_size(); i++) { + matrix_.increment_row_index(i); + matrix_.increment_col_index(i); + } + + // irn/jcn must be repointed even though structure is unchanged, since + // extractUpperTriangle (if symmetric) allocates a new matrix_ object + // with a new underlying buffer each call. + mumps_solver_.job = JOB_FACTORIZATION_PHASE; // reuse stored analysis/ordering + mumps_solver_.irn = matrix_.row_indices_data(); + mumps_solver_.jcn = matrix_.column_indices_data(); + mumps_solver_.a = matrix_.values_data(); + + dmumps_c(&mumps_solver_); + + if (INFOG(1) != 0) { + std::cerr << "MUMPS reported an error during factorization update " + << "(INFOG(1) = " << INFOG(1) << ").\n"; + } + + if (mumps_solver_.sym == SYM_POSITIVE_DEFINITE && INFOG(12) != 0) { + std::cerr << "Matrix declared positive definite, " + << "but negative pivots were encountered during refactorization " + << "(INFOG(12) = " << INFOG(12) << ").\n"; + } + } + private: void initialize(); void finalize(); diff --git a/include/LinearAlgebra/Solvers/csr_lu_solver.h b/include/LinearAlgebra/Solvers/csr_lu_solver.h index c5d32838e..1e6639cd8 100644 --- a/include/LinearAlgebra/Solvers/csr_lu_solver.h +++ b/include/LinearAlgebra/Solvers/csr_lu_solver.h @@ -85,6 +85,19 @@ class SparseLUSolver */ void solveInPlace(Vector& b) const; + /** + * @brief Recompute the LU factorization from a matrix with the same + * sparsity pattern as the one originally used to construct this solver. + * + * Reuses the existing RCM ordering (perm_/perm_inv_), since ordering + * depends only on the sparsity pattern, which is assumed unchanged here. + * Only the numeric factorization is redone. + * + * @param A Matrix with identical structure to the one originally + * factorized, but updated numeric values. + */ + void updateValues(const SparseMatrixCSR& A); + private: // LU decomposition data structures AllocatableVector L_values_, U_values_; // Non-zero values for L and U @@ -152,6 +165,23 @@ SparseLUSolver::SparseLUSolver(const SparseMatrixCSR +void SparseLUSolver::updateValues(const SparseMatrixCSR& A) +{ + assert(factorized_); + assert(A.rows() == A.columns()); + assert(A.rows() == static_cast(perm_.size()) && + "Matrix structure must match the originally factorized matrix"); + + // Reuse the existing RCM ordering - skips the O(n + nnz) RCM computation. + SparseMatrixCSR A_perm = permuteMatrix(A, perm_, perm_inv_); + + // Symbolic pattern (L_pattern/U_pattern) is identical to before since the + // sparsity pattern hasn't changed, so this recomputes the same patterns + // deterministically; only the numeric values differ downstream. + factorize(A_perm); +} + /** * Solves Ax = b for Vector type * @param b - Right-hand side vector (overwritten with solution) diff --git a/scripts/compile.sh b/scripts/compile.sh index b0537fb0a..5fe14b5d6 100755 --- a/scripts/compile.sh +++ b/scripts/compile.sh @@ -65,4 +65,4 @@ if [ -n "$build_type" ]; then .. || { echo "CMake configuration failed"; exit 1; } fi -cmake --build ${PWD}/../build -j 2 +cmake --build ${PWD}/../build -j 1 From 02e74b08167aa1033204e3a411d86f90870e605e Mon Sep 17 00:00:00 2001 From: julianlitz Date: Tue, 11 Aug 2026 22:58:56 +0200 Subject: [PATCH 2/2] correction --- include/DirectSolver/DirectSolverTake/directSolverTake.h | 5 +---- 1 file changed, 1 insertion(+), 4 deletions(-) diff --git a/include/DirectSolver/DirectSolverTake/directSolverTake.h b/include/DirectSolver/DirectSolverTake/directSolverTake.h index 4438f5268..603539a8c 100644 --- a/include/DirectSolver/DirectSolverTake/directSolverTake.h +++ b/include/DirectSolver/DirectSolverTake/directSolverTake.h @@ -22,11 +22,8 @@ class DirectSolverTake : public DirectSolver #else using SystemMatrix = SparseMatrixCSR; using SystemSolver = SparseLUSolver; - // 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: