Skip to content

Encapsulate MUMPS solver class #182

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

Open
julianlitz wants to merge 20 commits into
base: main
Choose a base branch
from
Open

Encapsulate MUMPS solver class #182

Show file tree
Hide file tree
Changes from all commits
Commits
Show all changes
20 commits
Select commit Hold shift + click to select a range
ff125ae
extract
julianlitz Feb 18, 2026
8bb865c
...
julianlitz Feb 18, 2026
83f0815
Reorder MUMPS solver structure declaration
julianlitz Feb 18, 2026
489375c
Reorganize MUMPS solver structure declaration
julianlitz Feb 18, 2026
ff96285
Implement CooMumpsSolver for sparse matrix solving
julianlitz Feb 19, 2026
f2b8165
Refactor CooMumpsSolver class interface
julianlitz Feb 19, 2026
a47f84b
Add linear algebra sources to CMakeLists
julianlitz Feb 19, 2026
88eb792
Expand COO_MUMPS_Solver
julianlitz Feb 21, 2026
697740a
Formatting
julianlitz Feb 21, 2026
7421a1e
Fix mumps solver test
julianlitz Feb 21, 2026
d93ccaf
Update coo_mumps_solver.cpp
julianlitz Feb 22, 2026
8f88e53
Update smootherGive.cpp
julianlitz Feb 23, 2026
6dcf5d0
Update smootherTake.cpp
julianlitz Feb 23, 2026
7f62372
Update extrapolatedSmootherGive.cpp
julianlitz Feb 23, 2026
fdc909f
Update extrapolatedSmootherTake.cpp
julianlitz Feb 23, 2026
7509f20
pragma once in .cpp tests file removed
julianlitz Mar 6, 2026
dd9eba4
avoid std::optional in DirectSolver
julianlitz Mar 6, 2026
ac69055
avoid std::optional in DirectSolver
julianlitz Mar 6, 2026
7582b61
Remove std::optional und separate inner matrix construction
julianlitz Mar 8, 2026
12c7c3f
Remove unused variables
julianlitz Mar 9, 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
6 changes: 3 additions & 3 deletions include/Definitions/geometry_helper.h
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -23,9 +23,9 @@ inline void compute_jacobian_elements(const DomainGeometry& domain_geometry, dou
/* which is represented by: */
/* [arr, 0.5*art] */
/* [0.5*atr, att] */
arr = 0.5 * (Jtt * Jtt + Jrt * Jrt) * coeff_alpha / fabs(detDF);
att = 0.5 * (Jtr * Jtr + Jrr * Jrr) * coeff_alpha / fabs(detDF);
art = (-Jtt * Jtr - Jrt * Jrr) * coeff_alpha / fabs(detDF);
arr = 0.5 * (Jtt * Jtt + Jrt * Jrt) * coeff_alpha / std::fabs(detDF);
att = 0.5 * (Jtr * Jtr + Jrr * Jrr) * coeff_alpha / std::fabs(detDF);
art = (-Jtt * Jtr - Jrt * Jrr) * coeff_alpha / std::fabs(detDF);
/* Note that the inverse Jacobian matrix DF^{-1} is: */
/* 1.0 / det(DF) * */
/* [Jtt, -Jrt] */
Expand Down
21 changes: 6 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 @@ -12,14 +12,12 @@ class DirectSolver_COO_MUMPS_Give : public DirectSolver
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 build in the member initializer of the DirectSolver class.

// clang-format off
const Stencil stencil_interior_ = {
Expand Down Expand Up @@ -49,15 +47,14 @@ class DirectSolver_COO_MUMPS_Give : public DirectSolver
};
// 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 @@ -69,12 +66,6 @@ class DirectSolver_COO_MUMPS_Give : public DirectSolver
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 @@ -89,4 +80,4 @@ class DirectSolver_COO_MUMPS_Give : public DirectSolver
double detDF, double coeff_beta);
};

#endif
#endif
21 changes: 6 additions & 15 deletions include/DirectSolver/DirectSolver-COO-MUMPS-Take/directSolverTake.h
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -12,14 +12,12 @@ class DirectSolver_COO_MUMPS_Take : public DirectSolver
const DensityProfileCoefficients& density_profile_coefficients,
bool DirBC_Interior, int num_omp_threads);

~DirectSolver_COO_MUMPS_Take() 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 build in the member initializer of the DirectSolver class.

// clang-format off
const Stencil stencil_interior_ = {
Expand Down Expand Up @@ -49,15 +47,14 @@ class DirectSolver_COO_MUMPS_Take : public DirectSolver
};
// 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 @@ -69,12 +66,6 @@ class DirectSolver_COO_MUMPS_Take : public DirectSolver
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 @@ -90,4 +81,4 @@ class DirectSolver_COO_MUMPS_Take : public DirectSolver
ConstVector<double>& coeff_beta);
};

#endif
#endif
1 change: 1 addition & 0 deletions include/DirectSolver/directSolver.h
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -17,6 +17,7 @@ class Level;
#include "../LinearAlgebra/Matrix/coo_matrix.h"
#include "../LinearAlgebra/Matrix/csr_matrix.h"
#include "../LinearAlgebra/Solvers/csr_lu_solver.h"
#include "../LinearAlgebra/Solvers/coo_mumps_solver.h"
#include "../Stencil/stencil.h"

#ifdef GMGPOLAR_USE_MUMPS
Expand Down
123 changes: 64 additions & 59 deletions include/ExtrapolatedSmoother/ExtrapolatedSmootherGive/extrapolatedSmootherGive.h
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -61,9 +61,6 @@ class ExtrapolatedSmootherGive : public ExtrapolatedSmoother
const DensityProfileCoefficients& density_profile_coefficients,
bool DirBC_Interior, int num_omp_threads);

// If MUMPS is enabled, this cleans up the inner boundary solver.
~ExtrapolatedSmootherGive() override;

// Performs one full coupled extrapolated smoothing sweep:
// BC -> WC -> BR -> WR
// Parallel implementation using OpenMP:
Expand All @@ -72,40 +69,13 @@ class ExtrapolatedSmootherGive : public ExtrapolatedSmoother
void extrapolatedSmoothing(Vector<double> x, ConstVector<double> rhs, Vector<double> temp) override;

private:
/* ------------------- */
/* Tridiagonal solvers */
/* ------------------- */

// Batched solver for cyclic-tridiagonal circle line A_sc matrices.
BatchedTridiagonalSolver<double> circle_tridiagonal_solver_;

// Batched solver for tridiagonal radial circle line A_sc matrices.
BatchedTridiagonalSolver<double> radial_tridiagonal_solver_;

// The A_sc matrix on i_r = 0 (inner circle) is NOT tridiagonal because
// it potentially includes across-origin coupling. Therefore, it is assembled
// into a sparse matrix and solved using a general-purpose sparse solver.
// When using the MUMPS solver, the matrix is assembled in COO format.
// When using the in-house solver, the matrix is stored in CSR format.
#ifdef GMGPOLAR_USE_MUMPS
using MatrixType = SparseMatrixCOO<double>;
DMUMPS_STRUC_C inner_boundary_mumps_solver_;
#else
using MatrixType = SparseMatrixCSR<double>;
SparseLUSolver<double> inner_boundary_lu_solver_;
#endif
// Sparse matrix for the non-tridiagonal inner boundary circle block.
MatrixType inner_boundary_circle_matrix_;

// Note:
// - circle_tridiagonal_solver_[batch=0] is unused. Use the COO/CSR matrix instead.
// - circle_tridiagonal_solver_[batch=i_r] solves circle line i_r.
// - radial_tridiagonal_solver_[batch=i_theta] solves radial line i_theta.

/* ------------------- */
/* Stencil definitions */
/* ------------------- */

// The stencil definitions must be defined before the declaration of the inner_boundary_mumps_solver_,
// since the mumps solver will be build in the member initializer of the Smoother class.

// Stencils encode neighborhood connectivity for A_sc matrix assembly.
// It is only used in the construction of COO/CSR matrices.
// Thus it is only used for the interior boundary matrix and not needed for the tridiagonal matrices.
Expand All @@ -127,6 +97,45 @@ class ExtrapolatedSmootherGive : public ExtrapolatedSmoother
};
// clang-format on

/* ------------------- */
/* Tridiagonal solvers */
/* ------------------- */

// Batched solver for cyclic-tridiagonal circle line A_sc matrices.
BatchedTridiagonalSolver<double> circle_tridiagonal_solver_;

// Batched solver for tridiagonal radial line A_sc matrices.
BatchedTridiagonalSolver<double> radial_tridiagonal_solver_;

// Note:
// - circle_tridiagonal_solver_[batch=0] is unused. Use the COO/CSR matrix instead.
// - circle_tridiagonal_solver_[batch=i_r] solves circle line i_r.
// - radial_tridiagonal_solver_[batch=i_theta] solves radial line i_theta.

/* ------------------------ */
/* Interior boundary solver */
/* ------------------------ */

// The A_sc matrix on i_r = 0 (inner circle) is NOT tridiagonal because
// it potentially includes across-origin coupling. Therefore, it is assembled
// into a sparse matrix and solved using a general-purpose sparse solver.
// When using the MUMPS solver, the matrix is assembled in COO format.
// When using the in-house solver, the matrix is stored in CSR format.

#ifdef GMGPOLAR_USE_MUMPS
// When using the MUMPS solver, the matrix is assembled in COO format.
using MatrixType = SparseMatrixCOO<double>;
// MUMPS solver structure with the solver matrix initialized in the constructor.
CooMumpsSolver inner_boundary_mumps_solver_;
#else
// When using the in-house solver, the matrix is stored in CSR format.
using MatrixType = SparseMatrixCSR<double>;
// Sparse matrix for the non-tridiagonal inner boundary circle block.
MatrixType inner_boundary_circle_matrix_;
// LU solver for the inner boundary circle block.
SparseLUSolver<double> inner_boundary_lu_solver_;
#endif

// Select correct stencil depending on the grid position.
const Stencil& getStencil(int i_r, int i_theta) const; /* Only i_r = 0 implemented */
// Number of nonzero A_sc entries.
Expand All @@ -138,30 +147,36 @@ class ExtrapolatedSmootherGive : public ExtrapolatedSmoother
/* --------------- */
/* Matrix assembly */
/* --------------- */

// Build all A_sc matrices for circle and radial smoothers.
void buildAscMatrices();
// Build A_sc matrix block for a single circular line.
void buildAscCircleSection(int i_r);
// Build A_sc matrix block for a single radial line.
void buildAscRadialSection(int i_theta);
// Build A_sc for a specific node (i_r, i_theta)
void nodeBuildAscGive(int i_r, int i_theta, const PolarGrid& grid, bool DirBC_Interior,
MatrixType& inner_boundary_circle_matrix,
BatchedTridiagonalSolver<double>& circle_tridiagonal_solver,
BatchedTridiagonalSolver<double>& radial_tridiagonal_solver, double arr, double att,
double art, double detDF, double coeff_beta);
void buildTridiagonalSolverMatrices();
void buildTridiagonalCircleSection(int i_r);
void buildTridiagonalRadialSection(int i_theta);
// Build the tridiagonal solver matrices for a specific node (i_r, i_theta)
void nodeBuildTridiagonalSolverMatrices(int i_r, int i_theta, const PolarGrid& grid, bool DirBC_Interior,
BatchedTridiagonalSolver<double>& circle_tridiagonal_solver,
BatchedTridiagonalSolver<double>& radial_tridiagonal_solver, double arr,
double att, double art, double detDF, double coeff_beta);

// Build the solver matrix for the interior boundary (i_r = 0) which is non-tridiagonal due to across-origin coupling.
MatrixType buildInteriorBoundarySolverMatrix();
// Build the solver matrix for a specific node (i_r = 0, i_theta) on the interior boundary.
void nodeBuildInteriorBoundarySolverMatrix_i_r_0(int i_theta, const PolarGrid& grid, bool DirBC_Interior,
MatrixType& matrix, double arr, double att, double art,
double detDF, double coeff_beta);
void nodeBuildInteriorBoundarySolverMatrix_i_r_1(int i_theta, const PolarGrid& grid, bool DirBC_Interior,
MatrixType& matrix, double arr, double att, double art,
double detDF, double coeff_beta);

/* ---------------------- */
/* Orthogonal application */
/* ---------------------- */

// Compute temp = f_sc − A_sc^ortho * u_sc^ortho (precomputed right-hand side)
// where x = u_sc and rhs = f_sc
void applyAscOrthoCircleSection(int i_r, SmootherColor smoother_color, ConstVector<double> x,
ConstVector<double> rhs, Vector<double> temp);
void applyAscOrthoRadialSection(int i_theta, SmootherColor smoother_color, ConstVector<double> x,
ConstVector<double> rhs, Vector<double> temp);
void applyAscOrthoBlackCircleSection(ConstVector<double> x, ConstVector<double> rhs, Vector<double> temp);
void applyAscOrthoWhiteCircleSection(ConstVector<double> x, ConstVector<double> rhs, Vector<double> temp);
void applyAscOrthoBlackRadialSection(ConstVector<double> x, ConstVector<double> rhs, Vector<double> temp);
void applyAscOrthoWhiteRadialSection(ConstVector<double> x, ConstVector<double> rhs, Vector<double> temp);

/* ----------------- */
/* Line-wise solvers */
Expand All @@ -179,14 +194,4 @@ class ExtrapolatedSmootherGive : public ExtrapolatedSmoother
void solveWhiteCircleSection(Vector<double> x, Vector<double> temp);
void solveBlackRadialSection(Vector<double> x, Vector<double> temp);
void solveWhiteRadialSection(Vector<double> x, Vector<double> temp);

/* ----------------------------------- */
/* Initialize and destroy MUMPS solver */
/* ----------------------------------- */
#ifdef GMGPOLAR_USE_MUMPS
// Initialize sparse MUMPS solver with assembled COO matrix.
void initializeMumpsSolver(DMUMPS_STRUC_C& mumps_solver, SparseMatrixCOO<double>& solver_matrix);
// Release MUMPS internal memory and MPI structures.
void finalizeMumpsSolver(DMUMPS_STRUC_C& mumps_solver);
#endif
};
Loading