Skip to content

Commit

Permalink
Support matrix-free for eigenvalue solvers from SLEPc
Browse files Browse the repository at this point in the history 
Closes idaholab#7398
  • Loading branch information
@fdkong
fdkong committed May 18, 2020
1 parent 7a584ba commit ed8f2ab
Show file tree
Hide file tree
Showing 7 changed files with 192 additions and 67 deletions.
7 changes: 7 additions & 0 deletions framework/include/problems/EigenProblem.h
View file
Original file line number Diff line number Diff line change
Expand Up @@ -32,6 +32,8 @@ class EigenProblem : public FEProblemBase

virtual void solve() override;

virtual void init() override;

virtual bool converged() override;

virtual unsigned int getNEigenPairsRequired() { return _n_eigen_pairs_required; }
Expand All @@ -55,6 +57,10 @@ class EigenProblem : public FEProblemBase
*/
bool negativeSignEigenKernel() { return _negative_sign_eigen_kernel; }

bool matrixFree() { return _matrix_free; }

void matrixFree(bool matrix_free) { _matrix_free = matrix_free; }

#if LIBMESH_HAVE_SLEPC
void setEigenproblemType(Moose::EigenProblemType eigen_problem_type);

Expand Down Expand Up @@ -103,6 +109,7 @@ class EigenProblem : public FEProblemBase
std::shared_ptr<NonlinearEigenSystem> _nl_eigen;

bool _negative_sign_eigen_kernel;
bool _matrix_free;

unsigned int _active_eigen_index;
/// Timers
Expand Down
2 changes: 2 additions & 0 deletions framework/include/systems/NonlinearEigenSystem.h
View file
Original file line number Diff line number Diff line change
Expand Up @@ -72,6 +72,8 @@ class NonlinearEigenSystem : public NonlinearSystemBase

virtual void initialSetup() override;

void attachSLEPcCallbacks();

/**
* Get the number of converged eigenvalues
*
Expand Down
7 changes: 7 additions & 0 deletions framework/include/utils/SlepcSupport.h
View file
Original file line number Diff line number Diff line change
Expand Up @@ -43,6 +43,13 @@ PetscErrorCode mooseSlepcEigenFormJacobianB(SNES snes, Vec x, Mat jac, Mat pc, v
PetscErrorCode mooseSlepcEigenFormFunctionA(SNES snes, Vec x, Vec r, void * ctx);
PetscErrorCode mooseSlepcEigenFormFunctionB(SNES snes, Vec x, Vec r, void * ctx);
PetscErrorCode mooseSlepcEigenFormFunctionAB(SNES snes, Vec x, Vec Ax, Vec Bx, void * ctx);

void attachCallbacksToMat(EigenProblem & eigen_problem, Mat mat, bool eigen);

PetscErrorCode mooseMatMult_Eigen(Mat mat, Vec x, Vec y);
PetscErrorCode mooseMatMult_NonEigen(Mat mat, Vec x, Vec y);
void setOperationsForShellMat(EigenProblem & eigen_problem, Mat mat, bool eigen);

} // namespace SlepcSupport
} // namespace moose

Expand Down
9 changes: 9 additions & 0 deletions framework/src/executioners/Eigenvalue.C
View file
Original file line number Diff line number Diff line change
Expand Up @@ -26,6 +26,13 @@ Eigenvalue::validParams()

params.addClassDescription("Eigenvalue solves a standard/generalized eigenvaue problem");

params.addParam<bool>(
"matrix_free",
false,
"Whether or not to use a matrix free fashion to form operators. "
"If true, shell matrices will be used and meanwhile a preconditioning matrix"
"may be formed as well.");

params.addPrivateParam<bool>("_use_eigen_value", true);

params.addParam<PostprocessorName>(
Expand Down Expand Up @@ -55,6 +62,8 @@ Eigenvalue::Eigenvalue(const InputParameters & parameters)

Moose::SlepcSupport::storeSlepcEigenProblemOptions(_eigen_problem, parameters);
_eigen_problem.setEigenproblemType(_eigen_problem.solverParams()._eigen_problem_type);

_eigen_problem.matrixFree(getParam<bool>("matrix_free"));
#endif

if (!parameters.isParamValid("normalization") && parameters.isParamSetByUser("normal_factor"))
Expand Down
16 changes: 14 additions & 2 deletions framework/src/problems/EigenProblem.C
View file
Original file line number Diff line number Diff line change
Expand Up @@ -289,8 +289,8 @@ EigenProblem::solve()
{
TIME_SECTION(_solve_timer);

_nl->solve();
_nl->update();
_nl_eigen->solve();
_nl_eigen->update();
}

// sync solutions in displaced problem
Expand All @@ -303,6 +303,18 @@ EigenProblem::solve()
#endif
}

void
EigenProblem::init()
{
// If matrix_free=true, this set Libmesh to use shell matrices
_nl_eigen->sys().use_shell_matrices(_matrix_free);

FEProblemBase::init();

// Set proper operations
_nl_eigen->attachSLEPcCallbacks();
}

bool
EigenProblem::converged()
{
Expand Down
122 changes: 61 additions & 61 deletions framework/src/systems/NonlinearEigenSystem.C
View file
Original file line number Diff line number Diff line change
Expand Up @@ -25,6 +25,7 @@
#include "libmesh/libmesh_config.h"
#include "libmesh/petsc_matrix.h"
#include "libmesh/sparse_matrix.h"
#include "libmesh/petsc_shell_matrix.h"

#if LIBMESH_HAVE_SLEPC

Expand All @@ -38,6 +39,15 @@ assemble_matrix(EquationSystems & es, const std::string & system_name)
EigenSystem & eigen_system = es.get_system<EigenSystem>(system_name);
NonlinearEigenSystem & eigen_nl = p->getNonlinearEigenSystem();

// If we use shell matrices, we only need to form a preconditioning matrix
if (eigen_system.use_shell_matrices())
{
p->computeJacobianTag(*eigen_system.current_local_solution.get(),
*eigen_system.precond_matrix,
eigen_nl.nonEigenMatrixTag());
return;
}

// If it is a linear generalized eigenvalue problem,
// we assemble A and B together
if (!p->isNonlinearEigenvalueSolver() && eigen_system.generalized())
Expand All @@ -54,72 +64,14 @@ assemble_matrix(EquationSystems & es, const std::string & system_name)
return;
}

// If it is a linear eigenvalue problem, we assemble matrix A
if (!p->isNonlinearEigenvalueSolver())
{
p->computeJacobianTag(*eigen_system.current_local_solution.get(),
*eigen_system.matrix_A,
eigen_nl.nonEigenMatrixTag());
}
else
{
Mat petsc_mat_A = static_cast<PetscMatrix<Number> &>(*eigen_system.matrix_A).mat();

PetscObjectComposeFunction((PetscObject)petsc_mat_A,
"formJacobian",
Moose::SlepcSupport::mooseSlepcEigenFormJacobianA);
PetscObjectComposeFunction((PetscObject)petsc_mat_A,
"formFunction",
Moose::SlepcSupport::mooseSlepcEigenFormFunctionA);

PetscObjectComposeFunction((PetscObject)petsc_mat_A,
"formFunctionAB",
Moose::SlepcSupport::mooseSlepcEigenFormFunctionAB);

PetscContainer container;
PetscContainerCreate(eigen_system.comm().get(), &container);
PetscContainerSetPointer(container, p);
PetscObjectCompose((PetscObject)petsc_mat_A, "formJacobianCtx", nullptr);
PetscObjectCompose((PetscObject)petsc_mat_A, "formJacobianCtx", (PetscObject)container);
PetscObjectCompose((PetscObject)petsc_mat_A, "formFunctionCtx", nullptr);
PetscObjectCompose((PetscObject)petsc_mat_A, "formFunctionCtx", (PetscObject)container);
PetscContainerDestroy(&container);

// Let libmesh do not close matrices before solve
eigen_system.eigen_solver->set_close_matrix_before_solve(false);
}
if (eigen_system.generalized())
{
if (eigen_system.matrix_B)
{
if (!p->isNonlinearEigenvalueSolver())
{
p->computeJacobianTag(*eigen_system.current_local_solution.get(),
*eigen_system.matrix_B,
eigen_nl.eigenMatrixTag());
}
else
{
Mat petsc_mat_B = static_cast<PetscMatrix<Number> &>(*eigen_system.matrix_B).mat();

PetscObjectComposeFunction((PetscObject)petsc_mat_B,
"formJacobian",
Moose::SlepcSupport::mooseSlepcEigenFormJacobianB);
PetscObjectComposeFunction((PetscObject)petsc_mat_B,
"formFunction",
Moose::SlepcSupport::mooseSlepcEigenFormFunctionB);

PetscContainer container;
PetscContainerCreate(eigen_system.comm().get(), &container);
PetscContainerSetPointer(container, p);
PetscObjectCompose((PetscObject)petsc_mat_B, "formFunctionCtx", nullptr);
PetscObjectCompose((PetscObject)petsc_mat_B, "formFunctionCtx", (PetscObject)container);
PetscObjectCompose((PetscObject)petsc_mat_B, "formJacobianCtx", nullptr);
PetscObjectCompose((PetscObject)petsc_mat_B, "formJacobianCtx", (PetscObject)container);
PetscContainerDestroy(&container);
}
}
else
mooseError("It is a generalized eigenvalue problem but matrix B is empty\n");
return;
}
}
}
Expand Down Expand Up @@ -198,7 +150,7 @@ NonlinearEigenSystem::solve()
// In DEBUG mode, Libmesh will check the residual automatically. This may cause
// an error because B does not need to assembly by default.
#ifdef DEBUG
if (_eigen_problem.isGeneralizedEigenvalueProblem())
if (_eigen_problem.isGeneralizedEigenvalueProblem() && sys().matrix_B)
sys().matrix_B->close();
#endif
// Solve the transient problem if we have a time integrator; the
Expand All @@ -224,6 +176,54 @@ NonlinearEigenSystem::solve()
_eigen_values[n] = getNthConvergedEigenvalue(n);
}

void
NonlinearEigenSystem::attachSLEPcCallbacks()
{
// Matrix A
if (_transient_sys.matrix_A)
{
Mat mat = static_cast<PetscMatrix<Number> &>(*_transient_sys.matrix_A).mat();

Moose::SlepcSupport::attachCallbacksToMat(_eigen_problem, mat, false);

// Let libmesh do not close matrices before solve
_transient_sys.eigen_solver->set_close_matrix_before_solve(false);
}

// Matrix B
if (_transient_sys.matrix_B)
{
Mat mat = static_cast<PetscMatrix<Number> &>(*_transient_sys.matrix_B).mat();

Moose::SlepcSupport::attachCallbacksToMat(_eigen_problem, mat, true);
}

// Shell matrix A
if (_transient_sys.shell_matrix_A)
{
Mat mat = static_cast<PetscShellMatrix<Number> &>(*_transient_sys.shell_matrix_A).mat();

// Attach callbacks for nonlinear eigenvalue solver
Moose::SlepcSupport::attachCallbacksToMat(_eigen_problem, mat, false);

// Set MatMult operations for shell
Moose::SlepcSupport::setOperationsForShellMat(_eigen_problem, mat, false);

_transient_sys.eigen_solver->set_close_matrix_before_solve(false);
}

// Shell matrix B
if (_transient_sys.shell_matrix_B)
{
Mat mat = static_cast<PetscShellMatrix<Number> &>(*_transient_sys.shell_matrix_B).mat();

Moose::SlepcSupport::attachCallbacksToMat(_eigen_problem, mat, true);

// Set MatMult operations for shell
Moose::SlepcSupport::setOperationsForShellMat(_eigen_problem, mat, true);
}
}

void
NonlinearEigenSystem::stopSolve()
{
Expand Down
96 changes: 92 additions & 4 deletions framework/src/utils/SlepcSupport.C
View file
Original file line number Diff line number Diff line change
Expand Up @@ -367,7 +367,6 @@ setEigenSolverOptions(SolverParams & solver_params, const InputParameters & para
Moose::PetscSupport::setSinglePetscOption("-eps_type", "power");
Moose::PetscSupport::setSinglePetscOption("-eps_power_nonlinear", "1");
Moose::PetscSupport::setSinglePetscOption("-eps_target_magnitude", "");
Moose::PetscSupport::setSinglePetscOption("-st_type", "sinvert");
#else
mooseError("Nonlinear Inverse Power requires SLEPc 3.7.3 or higher");
#endif
Expand All @@ -379,7 +378,6 @@ setEigenSolverOptions(SolverParams & solver_params, const InputParameters & para
Moose::PetscSupport::setSinglePetscOption("-eps_power_nonlinear", "1");
Moose::PetscSupport::setSinglePetscOption("-eps_power_snes_mf_operator", "1");
Moose::PetscSupport::setSinglePetscOption("-eps_target_magnitude", "");
Moose::PetscSupport::setSinglePetscOption("-st_type", "sinvert");
#else
mooseError("Matrix-free nonlinear Inverse Power requires SLEPc 3.7.3 or higher");
#endif
Expand All @@ -395,7 +393,6 @@ setEigenSolverOptions(SolverParams & solver_params, const InputParameters & para
Moose::PetscSupport::setSinglePetscOption(
"-init_eps_max_it", stringify(params.get<unsigned int>("free_power_iterations")));
Moose::PetscSupport::setSinglePetscOption("-eps_target_magnitude", "");
Moose::PetscSupport::setSinglePetscOption("-st_type", "sinvert");
#else
mooseError("Newton-based eigenvalue solver requires SLEPc 3.7.3 or higher");
#endif
Expand All @@ -413,7 +410,6 @@ setEigenSolverOptions(SolverParams & solver_params, const InputParameters & para
Moose::PetscSupport::setSinglePetscOption(
"-init_eps_max_it", stringify(params.get<unsigned int>("free_power_iterations")));
Moose::PetscSupport::setSinglePetscOption("-eps_target_magnitude", "");
Moose::PetscSupport::setSinglePetscOption("-st_type", "sinvert");
#else
mooseError("Matrix-free Newton-based eigenvalue solver requires SLEPc 3.7.3 or higher");
#endif
Expand Down Expand Up @@ -568,6 +564,98 @@ mooseSlepcEigenFormFunctionAB(SNES /*snes*/, Vec x, Vec Ax, Vec Bx, void * ctx)
PetscFunctionReturn(0);
}

void
attachCallbacksToMat(EigenProblem & eigen_problem, Mat mat, bool eigen)
{
PetscObjectComposeFunction((PetscObject)mat,
"formJacobian",
eigen ? Moose::SlepcSupport::mooseSlepcEigenFormJacobianB
: Moose::SlepcSupport::mooseSlepcEigenFormJacobianA);
PetscObjectComposeFunction((PetscObject)mat,
"formFunction",
eigen ? Moose::SlepcSupport::mooseSlepcEigenFormFunctionB
: Moose::SlepcSupport::mooseSlepcEigenFormFunctionA);

PetscObjectComposeFunction(
(PetscObject)mat, "formFunctionAB", Moose::SlepcSupport::mooseSlepcEigenFormFunctionAB);

PetscContainer container;
PetscContainerCreate(eigen_problem.comm().get(), &container);
PetscContainerSetPointer(container, &eigen_problem);
PetscObjectCompose((PetscObject)mat, "formJacobianCtx", nullptr);
PetscObjectCompose((PetscObject)mat, "formJacobianCtx", (PetscObject)container);
PetscObjectCompose((PetscObject)mat, "formFunctionCtx", nullptr);
PetscObjectCompose((PetscObject)mat, "formFunctionCtx", (PetscObject)container);
PetscContainerDestroy(&container);
}

void
mooseMatMult(EigenProblem & eigen_problem, Vec x, Vec r, TagID tag)
{
NonlinearSystemBase & nl = eigen_problem.getNonlinearSystemBase();
System & sys = nl.system();

PetscVector<Number> X_global(x, sys.comm()), R(r, sys.comm());

// update local solution
X_global.localize(*sys.current_local_solution.get());

R.zero();

eigen_problem.computeResidualTag(*sys.current_local_solution.get(), R, tag);

R.close();
}

PetscErrorCode
mooseMatMult_Eigen(Mat mat, Vec x, Vec r)
{
PetscFunctionBegin;
void * ctx = nullptr;
MatShellGetContext(mat, &ctx);

if (!ctx)
mooseError("No context is set for shell matrix ");

EigenProblem * eigen_problem = static_cast<EigenProblem *>(ctx);
NonlinearEigenSystem & eigen_nl = eigen_problem->getNonlinearEigenSystem();

mooseMatMult(*eigen_problem, x, r, eigen_nl.eigenVectorTag());

if (eigen_problem->negativeSignEigenKernel())
VecScale(r, -1.);

PetscFunctionReturn(0);
}

PetscErrorCode
mooseMatMult_NonEigen(Mat mat, Vec x, Vec r)
{
PetscFunctionBegin;
void * ctx = nullptr;
MatShellGetContext(mat, &ctx);

if (!ctx)
mooseError("No context is set for shell matrix ");

EigenProblem * eigen_problem = static_cast<EigenProblem *>(ctx);
NonlinearEigenSystem & eigen_nl = eigen_problem->getNonlinearEigenSystem();

mooseMatMult(*eigen_problem, x, r, eigen_nl.nonEigenVectorTag());

PetscFunctionReturn(0);
}

void
setOperationsForShellMat(EigenProblem & eigen_problem, Mat mat, bool eigen)
{
MatShellSetContext(mat, &eigen_problem);
MatShellSetOperation(mat,
MATOP_MULT,
eigen ? (void (*)(void))mooseMatMult_Eigen
: (void (*)(void))mooseMatMult_NonEigen);
}

} // namespace SlepcSupport
} // namespace moose

Expand Down

0 comments on commit ed8f2ab

Please sign in to comment.