[LinearSolvers] Adding CG factorization for Spectra #13619

takeshimg92 · 2025-07-11T04:08:53Z

📝 Description

OBS: not sure if related to #13587

Problem

Kratos crashes and returns a "Killed" message when trying to calculate the eigenvalues for structural mechanics problems for very large meshes.

Cause

By default, Kratosuses direct solvers to perform matrix operations in eigenvalue problems.

Fix

Modified spectra_sym_g_eigs_shift_solver.h in LinearSolversApplication to make it work with large meshes. Now the user can add a custom setting in the ProjectParameters.json file when using the "solver_type": "spectra_sym_g_eigs_shift". Example:

        "eigensolver_settings": {
            "solver_type": "spectra_sym_g_eigs_shift",
            "max_iteration": 500,
            "shift": 0.0,
            "number_of_eigenvalues": 5,
            "echo_level": 1,
            "custom_factorization": "sparse_lu"
       },

"custom_factorization": "sparse_lu" was the (implicit) default choice implemented in Kratos. This employs the sparse_lu direct solver to perform operations with matrices. However, direct methods require a lot of memory, causing the simulations to get killed when using large meshes (even with only $\mathcal{O}(10^5)$ nodes). However, now the user can set "custom_factorization": "CG", which allows to use an iterative solver to perform such operations, reducing significantly the amount of memory required. Example:

        "eigensolver_settings": {
            "solver_type": "spectra_sym_g_eigs_shift",
            "max_iteration": 500,
            "shift": 0.0,
            "number_of_eigenvalues": 5,
            "echo_level": 1,
            "custom_factorization": "CG"
       },

For now these are the only two options.

Recommendations

For small meshes ($\lesssim \mathcal{O}(10^5)$ nodes), the option "custom_factorization": "sparse_lu" is the quickest; for large meshes ($\mathcal{O}(10^5)$ nodes) and above, use "custom_factorization": "CG".

roigcarlo

Probably not the best person to review this but just a couple of comments.

If we merge this, please, be sure to update these files with the new configuration option:

Kratos/applications/LinearSolversApplication/README.md
Kratos/applications/LinearSolversApplication/tests/test_eigensystem_solver.py
Kratos/applications/StructuralMechanicsApplication/python_scripts/structural_mechanics_eigensolver.py

Or make custom_factorization parameter optional. Otherwise it will crash on
https://github.com/KratosMultiphysics/Kratos/pull/13619/files#diff-a3d66cf319f1be32a63cf1686e3a86a71c328e4798458a27ba83550f55ae19a1R103.

On a side note #13587 is still not merged, so there should not be any effect from that PR.

roigcarlo · 2025-07-11T08:18:59Z

applications/LinearSolversApplication/custom_solvers/spectra_sym_g_eigs_shift_solver.h

+            if (m_cust_fact == "CG"){
+
+                MapConstVec x(x_in, m_n);
+                MapVec y(y_out, m_n);
+                y.noalias() = m_solverCG.solve(x);
+
+            } else {
+                MapConstVec x(x_in, m_n);
+                MapVec y(y_out, m_n);
+                y.noalias() = m_solver.solve(x);
+            }


x and y declarations can remain outside the branch

roigcarlo · 2025-07-11T08:19:05Z

applications/LinearSolversApplication/custom_solvers/spectra_sym_g_eigs_shift_solver.h


 namespace Kratos
 {

+class OwnSymShiftInvert;


Why do you need a forward declaration here?

matekelemen

I like the idea, but I think it would benefit from more generality than hard-coding CG.

I don't have much experience with eigensolvers so call me out if I'm writing something stupid. What I'd do is use our LinearSolver interface and its related factory

Kratos/kratos/solving_strategies/builder_and_solvers/p_multigrid/p_grid.cpp

Lines 87 to 93 in c69adba

    
           const std::string solver_name = SmootherSettings["solver_type"].Get<std::string>(); 
        
           using SolverFactoryRegistry = KratosComponents<LinearSolverFactory<TSparse,TDense>>; 
        
           KRATOS_ERROR_IF_NOT(SolverFactoryRegistry::Has(solver_name)) 
        
               << "\"" << solver_name << "\" is not a valid linear solver name in the registry. " 
        
               << "Make sure you imported the application it is defined in and that the spelling is correct."; 
        
           const auto& r_factory = SolverFactoryRegistry::Get(solver_name); 
        
           mpSolver = r_factory.Create(SmootherSettings);

LinearSolver::InitializeSolutionStep is supposed to do factorization-like tasks that can later be reused for different right hand sides with LinearSolver::PerformSolutionStep.

P.S.: the benefit is that you could use any linear solver we have in Kratos instead of just CG or LU.

matekelemen · 2025-07-11T07:49:08Z

applications/LinearSolversApplication/custom_solvers/spectra_sym_g_eigs_shift_solver.h

+#if !defined(BOOST_BIND_GLOBAL_PLACEHOLDERS)
+#define BOOST_BIND_GLOBAL_PLACEHOLDERS


What's this for?

Oh this is for suppressing some warnings in old versions of boosts (1.63 to 1.74 or so), but agree, that should be defined globally in any case, not per file.

matekelemen · 2025-07-11T07:50:23Z

applications/LinearSolversApplication/custom_solvers/spectra_sym_g_eigs_shift_solver.h

@@ -24,15 +27,22 @@
 #include "linear_solvers/iterative_solver.h"
 #include "custom_utilities/ublas_wrapper.h"
 #include "utilities/builtin_timer.h"
+#include "linear_solvers/amgcl_solver.h"


That's a chunky include, and as far as I can tell you don't need it here.

matekelemen · 2025-07-11T08:00:14Z

applications/LinearSolversApplication/custom_solvers/spectra_sym_g_eigs_shift_solver.h


 namespace Kratos
 {

+class OwnSymShiftInvert;


Where's this forward declaration needed?

philbucher · 2025-07-11T08:43:49Z

Cool stuff, I will comment probably on the weekend, have a bunch of comments and suggestions

philbucher

Global question aside from the one regarding the factory:
Afaik the advantage of a direct solver is also when computing many eigenvalues. The factorization can be done once and reused, whereas with an iterative solver you have to recompute the entire solution.

How many eigenvalues did you compute in your comparisons?

I have not yet done a detailed review because I would like to clarify the global questions first

philbucher · 2025-07-13T21:06:50Z

applications/LinearSolversApplication/custom_solvers/spectra_sym_g_eigs_shift_solver.h

@@ -60,7 +70,8 @@ class SpectraSymGEigsShiftSolver
            "normalize_eigenvectors": false,
            "max_iteration": 1000,
            "shift": 0.0,
-            "echo_level": 1
+            "echo_level": 1,
+            "custom_factorization": "sparse_lu"


why only limit this to two solvers? I.e. you can construct a linear solver via the factory, thus leaving it completely up to the user

Basically same comment as @matekelemen

philbucher · 2025-07-13T21:12:09Z

I like the idea, but I think it would benefit from more generality than hard-coding CG.

I don't have much experience with eigensolvers so call me out if I'm writing something stupid. What I'd do is use our LinearSolver interface and its related factory

Kratos/kratos/solving_strategies/builder_and_solvers/p_multigrid/p_grid.cpp

Lines 87 to 93 in c69adba

const std::string solver_name = SmootherSettings["solver_type"].Get<std::string>();

using SolverFactoryRegistry = KratosComponents<LinearSolverFactory<TSparse,TDense>>;

KRATOS_ERROR_IF_NOT(SolverFactoryRegistry::Has(solver_name))

<< "\"" << solver_name << "\" is not a valid linear solver name in the registry. "

<< "Make sure you imported the application it is defined in and that the spelling is correct.";

const auto& r_factory = SolverFactoryRegistry::Get(solver_name);

mpSolver = r_factory.Create(SmootherSettings);

LinearSolver::InitializeSolutionStep is supposed to do factorization-like tasks that can later be reused for different right hand sides with LinearSolver::PerformSolutionStep.

P.S.: the benefit is that you could use any linear solver we have in Kratos instead of just CG or LU.

Check out how this is done here:

Kratos/kratos/linear_solvers/monotonicity_preserving_solver.h

Lines 114 to 118 in a83ea2a

    
                   Parameters default_parameters( R"( 
        
                   { 
        
                       "solver_type"                    : "monotonicity_preserving", 
        
                       "inner_solver_settings"          : { 
        
                           "preconditioner_type"            : "amg",

matekelemen · 2025-07-13T21:45:49Z

Afaik the advantage of a direct solver is also when computing many eigenvalues.

I guess the problem here is not performance but feasibility. He wrote he ran out of memory, in which case an iterative solver really is the only way to go, even if it scales worse with the number of requested eigenvalues.

matekelemen · 2025-07-14T12:22:32Z

OBS: not sure if related to #13587

If you make the changes to support LinearSolver, you'll eventually be able to use those solvers as well.

That said, they're unlikely to help you much in this case. Sure they're pretty efficient for the methods they implement, but they still have the memory complexity inherent to direct methods.

For large systems, I recommend you give AMGCLSolver a try. Right now it does not store its "factorization" but #13624 addresses this issue.

takeshimg92 · 2025-07-15T01:06:02Z

Thank you for the answers.
@matekelemen and @philbucher thanks for the reference to LinearSolver interface and Kratos/kratos/linear_solvers/monotonicity_preserving_solver.h. Will take a look (and learn from it) and see whether I can adapt.

implemented CG factorization

835e256

takeshimg92 requested review from a team as code owners July 11, 2025 04:08

takeshimg92 changed the title ~~implemented CG factorization~~ [LinearSolvers] Adding CG factorization for Spectra Jul 11, 2025

roigcarlo reviewed Jul 11, 2025

View reviewed changes

matekelemen reviewed Jul 11, 2025

View reviewed changes

philbucher reviewed Jul 13, 2025

View reviewed changes

Provide feedback

Saved searches

Use saved searches to filter your results more quickly

Uh oh!

[LinearSolvers] Adding CG factorization for Spectra #13619

[LinearSolvers] Adding CG factorization for Spectra #13619

Uh oh!

takeshimg92 commented Jul 11, 2025 •

edited

Loading

Uh oh!

roigcarlo left a comment

Uh oh!

roigcarlo Jul 11, 2025 •

edited

Loading

Uh oh!

roigcarlo Jul 11, 2025 •

edited

Loading

Uh oh!

matekelemen left a comment •

edited

Loading

Uh oh!

matekelemen Jul 11, 2025

Uh oh!

roigcarlo Jul 11, 2025

Uh oh!

matekelemen Jul 11, 2025

Uh oh!

matekelemen Jul 11, 2025

Uh oh!

philbucher commented Jul 11, 2025

Uh oh!

philbucher left a comment

Uh oh!

philbucher Jul 13, 2025

Uh oh!

philbucher commented Jul 13, 2025

Uh oh!

matekelemen commented Jul 13, 2025

Uh oh!

matekelemen commented Jul 14, 2025

Uh oh!

takeshimg92 commented Jul 15, 2025

Uh oh!

Uh oh!

	const std::string solver_name = SmootherSettings["solver_type"].Get<std::string>();
	using SolverFactoryRegistry = KratosComponents<LinearSolverFactory<TSparse,TDense>>;
	KRATOS_ERROR_IF_NOT(SolverFactoryRegistry::Has(solver_name))
	<< "\"" << solver_name << "\" is not a valid linear solver name in the registry. "
	<< "Make sure you imported the application it is defined in and that the spelling is correct.";
	const auto& r_factory = SolverFactoryRegistry::Get(solver_name);
	mpSolver = r_factory.Create(SmootherSettings);

		#if !defined(BOOST_BIND_GLOBAL_PLACEHOLDERS)
		#define BOOST_BIND_GLOBAL_PLACEHOLDERS

[LinearSolvers] Adding CG factorization for Spectra #13619

Are you sure you want to change the base?

[LinearSolvers] Adding CG factorization for Spectra #13619

Uh oh!

Conversation

takeshimg92 commented Jul 11, 2025 • edited Loading Uh oh! There was an error while loading. Please reload this page.

Uh oh!

Problem

Cause

Fix

Recommendations

Uh oh!

roigcarlo left a comment

Choose a reason for hiding this comment

Uh oh!

roigcarlo Jul 11, 2025 • edited Loading Uh oh! There was an error while loading. Please reload this page.

Uh oh!

Choose a reason for hiding this comment

Uh oh!

roigcarlo Jul 11, 2025 • edited Loading Uh oh! There was an error while loading. Please reload this page.

Uh oh!

Choose a reason for hiding this comment

Uh oh!

matekelemen left a comment • edited Loading Uh oh! There was an error while loading. Please reload this page.

Uh oh!

Choose a reason for hiding this comment

Uh oh!

matekelemen Jul 11, 2025

Choose a reason for hiding this comment

Uh oh!

roigcarlo Jul 11, 2025

Choose a reason for hiding this comment

Uh oh!

matekelemen Jul 11, 2025

Choose a reason for hiding this comment

Uh oh!

matekelemen Jul 11, 2025

Choose a reason for hiding this comment

Uh oh!

philbucher commented Jul 11, 2025

Uh oh!

philbucher left a comment

Choose a reason for hiding this comment

Uh oh!

philbucher Jul 13, 2025

Choose a reason for hiding this comment

Uh oh!

philbucher commented Jul 13, 2025

Uh oh!

matekelemen commented Jul 13, 2025

Uh oh!

matekelemen commented Jul 14, 2025

Uh oh!

takeshimg92 commented Jul 15, 2025

Uh oh!

Uh oh!

takeshimg92 commented Jul 11, 2025 •

edited

Loading

roigcarlo Jul 11, 2025 •

edited

Loading

roigcarlo Jul 11, 2025 •

edited

Loading

matekelemen left a comment •

edited

Loading