Skip to content

Substitute COO_Mumps_Solver in DirectSolver #194

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 9 commits into from
Mar 18, 2026
Merged

Substitute COO_Mumps_Solver in DirectSolver #194

Show file tree
Hide file tree
Changes from all commits
Commits
Show all changes
9 commits
Select commit Hold shift + click to select a range
bff90fd
Add GG Anderson Mumps Solver
julianlitz Mar 13, 2026
b22e53e
GG Anderson Formatting
julianlitz Mar 13, 2026
1a7393e
GG Anderson Copy Copy Cat
julianlitz Mar 13, 2026
015a1f9
Unlimited power!
julianlitz Mar 14, 2026
b5c37ad
Yeet
julianlitz Mar 14, 2026
b797925
Resolve compile error
julianlitz Mar 14, 2026
931305a
Update buildSolverMatrix.inl
julianlitz Mar 14, 2026
742ce79
Grammar
julianlitz Mar 18, 2026
3ae9fb3
Grammar
julianlitz Mar 18, 2026
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
18 changes: 4 additions & 14 deletions include/DirectSolver/DirectSolver-COO-MUMPS-Give/applySymmetryShift.inl
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -131,25 +131,15 @@ void DirectSolver_COO_MUMPS_Give<DomainGeometry>::applySymmetryShift(Vector<doub
{
const PolarGrid& grid = DirectSolver<DomainGeometry>::grid_;
const bool DirBC_Interior = DirectSolver<DomainGeometry>::DirBC_Interior_;
const int num_omp_threads = DirectSolver<DomainGeometry>::num_omp_threads_;

assert(std::ssize(x) == grid.numberOfNodes());
assert(grid.nr() >= 4);

#pragma omp parallel sections num_threads(num_omp_threads)
{
#pragma omp section
{
if (DirBC_Interior) {
applySymmetryShiftInnerBoundary(x);
}
}

#pragma omp section
{
applySymmetryShiftOuterBoundary(x);
}
if (DirBC_Interior) {
applySymmetryShiftInnerBoundary(x);
}

applySymmetryShiftOuterBoundary(x);
}

#endif
66 changes: 15 additions & 51 deletions include/DirectSolver/DirectSolver-COO-MUMPS-Give/buildSolverMatrix.inl
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -785,6 +785,7 @@ void DirectSolver_COO_MUMPS_Give<DomainGeometry>::buildSolverMatrixCircleSection
{
const PolarGrid& grid = DirectSolver<DomainGeometry>::grid_;
const LevelCache<DomainGeometry>& level_cache = DirectSolver<DomainGeometry>::level_cache_;
const bool DirBC_Interior = DirectSolver<DomainGeometry>::DirBC_Interior_;

const double r = grid.radius(i_r);
for (int i_theta = 0; i_theta < grid.ntheta(); i_theta++) {
Expand All @@ -795,8 +796,7 @@ void DirectSolver_COO_MUMPS_Give<DomainGeometry>::buildSolverMatrixCircleSection
level_cache.obtainValues(i_r, i_theta, global_index, r, theta, coeff_beta, arr, att, art, detDF);

// Build solver matrix at the current node
nodeBuildSolverMatrixGive(i_r, i_theta, grid, DirectSolver<DomainGeometry>::DirBC_Interior_, solver_matrix, arr,
att, art, detDF, coeff_beta);
nodeBuildSolverMatrixGive(i_r, i_theta, grid, DirBC_Interior, solver_matrix, arr, att, art, detDF, coeff_beta);
}
}

Expand All @@ -806,6 +806,7 @@ void DirectSolver_COO_MUMPS_Give<DomainGeometry>::buildSolverMatrixRadialSection
{
const PolarGrid& grid = DirectSolver<DomainGeometry>::grid_;
const LevelCache<DomainGeometry>& level_cache = DirectSolver<DomainGeometry>::level_cache_;
const bool DirBC_Interior = DirectSolver<DomainGeometry>::DirBC_Interior_;

const double theta = grid.theta(i_theta);
for (int i_r = grid.numberSmootherCircles(); i_r < grid.nr(); i_r++) {
Expand All @@ -816,32 +817,21 @@ void DirectSolver_COO_MUMPS_Give<DomainGeometry>::buildSolverMatrixRadialSection
level_cache.obtainValues(i_r, i_theta, global_index, r, theta, coeff_beta, arr, att, art, detDF);

// Build solver matrix at the current node
nodeBuildSolverMatrixGive(i_r, i_theta, grid, DirectSolver<DomainGeometry>::DirBC_Interior_, solver_matrix, arr,
att, art, detDF, coeff_beta);
nodeBuildSolverMatrixGive(i_r, i_theta, grid, DirBC_Interior, solver_matrix, arr, att, art, detDF, coeff_beta);
}
}

// clang-format off

/* ------------------------------------------------------------------------ */
/* If the indexing is not smoother-based, please adjust the access patterns */
template <concepts::DomainGeometry DomainGeometry>
SparseMatrixCOO<double> DirectSolver_COO_MUMPS_Give<DomainGeometry>::buildSolverMatrix()
{
const PolarGrid& grid = DirectSolver<DomainGeometry>::grid_;
const int num_omp_threads = DirectSolver<DomainGeometry>::num_omp_threads_;
const PolarGrid& grid = DirectSolver<DomainGeometry>::grid_;
const int num_omp_threads = DirectSolver<DomainGeometry>::num_omp_threads_;

const int n = grid.numberOfNodes();
const int nnz = getNonZeroCountSolverMatrix();

// Although the matrix is symmetric, we need to store all its entries, so we disable the symmetry.
SparseMatrixCOO<double> solver_matrix(n, n, nnz);
solver_matrix.is_symmetric(false);

#pragma omp parallel for num_threads(num_omp_threads) if (nnz > 10'000)
for (int i = 0; i < nnz; i++) {
solver_matrix.value(i) = 0.0;
}
solver_matrix.is_symmetric(true);

if (num_omp_threads == 1) {
/* Single-threaded execution */
Expand All @@ -858,25 +848,25 @@ SparseMatrixCOO<double> DirectSolver_COO_MUMPS_Give<DomainGeometry>::buildSolver
const int additional_radial_tasks = grid.ntheta() % 3;
const int num_radial_tasks = grid.ntheta() - additional_radial_tasks;

#pragma omp parallel num_threads(num_omp_threads)
#pragma omp parallel num_threads(num_omp_threads)
{
#pragma omp for
#pragma omp for
for (int circle_task = 0; circle_task < num_circle_tasks; circle_task += 3) {
int i_r = grid.numberSmootherCircles() - circle_task - 1;
buildSolverMatrixCircleSection(i_r, solver_matrix);
}
#pragma omp for
#pragma omp for
for (int circle_task = 1; circle_task < num_circle_tasks; circle_task += 3) {
int i_r = grid.numberSmootherCircles() - circle_task - 1;
buildSolverMatrixCircleSection(i_r, solver_matrix);
}
#pragma omp for nowait
#pragma omp for nowait
for (int circle_task = 2; circle_task < num_circle_tasks; circle_task += 3) {
int i_r = grid.numberSmootherCircles() - circle_task - 1;
buildSolverMatrixCircleSection(i_r, solver_matrix);
}

#pragma omp for
#pragma omp for
for (int radial_task = 0; radial_task < num_radial_tasks; radial_task += 3) {
if (radial_task > 0) {
int i_theta = radial_task + additional_radial_tasks;
Expand All @@ -892,7 +882,7 @@ SparseMatrixCOO<double> DirectSolver_COO_MUMPS_Give<DomainGeometry>::buildSolver
}
}
}
#pragma omp for
#pragma omp for
for (int radial_task = 1; radial_task < num_radial_tasks; radial_task += 3) {
if (radial_task > 1) {
int i_theta = radial_task + additional_radial_tasks;
Expand All @@ -911,41 +901,15 @@ SparseMatrixCOO<double> DirectSolver_COO_MUMPS_Give<DomainGeometry>::buildSolver
}
}
}
#pragma omp for
#pragma omp for
for (int radial_task = 2; radial_task < num_radial_tasks; radial_task += 3) {
int i_theta = radial_task + additional_radial_tasks;
buildSolverMatrixRadialSection(i_theta, solver_matrix);
}
}
}

/* Mumps: In the case of symmetric matrices, only half of the matrix should be provided. */
/* Speeds up factorization time. */
const bool construct_symmetric = true;

if (!construct_symmetric) {
return solver_matrix;
}

/* Only store the upper tridiagonal entries of the symmetric solver_matrix */
const int symmetric_nnz = nnz - (nnz - n) / 2;
SparseMatrixCOO<double> symmetric_solver_matrix(n, n, symmetric_nnz);
symmetric_solver_matrix.is_symmetric(true);

int current_nz = 0;
for (int nz_index = 0; nz_index < nnz; nz_index++) {
const int row = solver_matrix.row_index(nz_index);
const int col = solver_matrix.col_index(nz_index);
if (row <= col) {
symmetric_solver_matrix.row_index(current_nz) = row;
symmetric_solver_matrix.col_index(current_nz) = col;
symmetric_solver_matrix.value(current_nz) = std::move(solver_matrix.value(nz_index));
current_nz++;
}
}

return symmetric_solver_matrix;
return solver_matrix;
}
// clang-format on

#endif
20 changes: 5 additions & 15 deletions include/DirectSolver/DirectSolver-COO-MUMPS-Give/directSolverGive.h
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -13,14 +13,12 @@ class DirectSolver_COO_MUMPS_Give : public DirectSolver<DomainGeometry>
const DensityProfileCoefficients& density_profile_coefficients,
bool DirBC_Interior, int num_omp_threads);

~DirectSolver_COO_MUMPS_Give() override;
// Note: The rhs (right-hand side) vector gets overwritten during the solution process.
void solveInPlace(Vector<double> solution) override;

private:
// Solver matrix and MUMPS solver structure
SparseMatrixCOO<double> solver_matrix_;
DMUMPS_STRUC_C mumps_solver_;
// The stencil definitions must be defined before the declaration of the mumps_solver,
// since the mumps solver will be built in the member initializer of the DirectSolver class.

// clang-format off
const Stencil stencil_interior_ = {
Expand Down Expand Up @@ -50,15 +48,14 @@ class DirectSolver_COO_MUMPS_Give : public DirectSolver<DomainGeometry>
};
// clang-format on

// MUMPS solver structure with the solver matrix initialized in the constructor.
CooMumpsSolver mumps_solver_;

// Constructs a symmetric solver matrix.
SparseMatrixCOO<double> buildSolverMatrix();
void buildSolverMatrixCircleSection(const int i_r, SparseMatrixCOO<double>& solver_matrix);
void buildSolverMatrixRadialSection(const int i_theta, SparseMatrixCOO<double>& solver_matrix);

// Initializes the MUMPS solver with the specified matrix.
// Converts to 1-based indexing.
void initializeMumpsSolver(DMUMPS_STRUC_C& mumps_solver, SparseMatrixCOO<double>& solver_matrix);

// Adjusts the right-hand side vector for symmetry corrections.
// This modifies the system from
// A * solution = rhs
Expand All @@ -70,12 +67,6 @@ class DirectSolver_COO_MUMPS_Give : public DirectSolver<DomainGeometry>
void applySymmetryShiftInnerBoundary(Vector<double> x) const;
void applySymmetryShiftOuterBoundary(Vector<double> x) const;

// Solves the adjusted system symmetric(matrixA) * solution = rhs using the MUMPS solver.
void solveWithMumps(Vector<double> solution);

// Finalizes the MUMPS solver, releasing any allocated resources.
void finalizeMumpsSolver(DMUMPS_STRUC_C& mumps_solver);

// Returns the total number of non-zero elements in the solver matrix.
int getNonZeroCountSolverMatrix() const;

Expand All @@ -93,7 +84,6 @@ class DirectSolver_COO_MUMPS_Give : public DirectSolver<DomainGeometry>
#include "applySymmetryShift.inl"
#include "buildSolverMatrix.inl"
#include "directSolverGive.inl"
#include "initializeMumps.inl"
#include "matrixStencil.inl"

#endif
11 changes: 2 additions & 9 deletions include/DirectSolver/DirectSolver-COO-MUMPS-Give/directSolverGive.inl
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -8,9 +8,8 @@ DirectSolver_COO_MUMPS_Give<DomainGeometry>::DirectSolver_COO_MUMPS_Give(
const DensityProfileCoefficients& density_profile_coefficients, bool DirBC_Interior, int num_omp_threads)
: DirectSolver<DomainGeometry>(grid, level_cache, domain_geometry, density_profile_coefficients, DirBC_Interior,
num_omp_threads)
, mumps_solver_(buildSolverMatrix())
{
solver_matrix_ = buildSolverMatrix();
initializeMumpsSolver(mumps_solver_, solver_matrix_);
}

template <concepts::DomainGeometry DomainGeometry>
Expand All @@ -23,13 +22,7 @@ void DirectSolver_COO_MUMPS_Give<DomainGeometry>::solveInPlace(Vector<double> so
// 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.
solveWithMumps(solution);
}

template <concepts::DomainGeometry DomainGeometry>
DirectSolver_COO_MUMPS_Give<DomainGeometry>::~DirectSolver_COO_MUMPS_Give()
{
finalizeMumpsSolver(mumps_solver_);
mumps_solver_.solve(solution);
}

#endif
128 changes: 0 additions & 128 deletions include/DirectSolver/DirectSolver-COO-MUMPS-Give/initializeMumps.inl
View file
Open in desktop

This file was deleted.

Loading
Loading