Detailed Description

A class of linear solver based on PETSc routines.

All command-line options that are not recognized by OGS are passed on to PETSc, i.e., they potentially affect the linear solver. The linear solver options in the project file take precedence over the command-line options, because the former are processed at a later stage.

Definition at line 25 of file PETScLinearSolver.h.

#include <PETScLinearSolver.h>

Public Member Functions
	PETScLinearSolver (std::string const &prefix, std::string const &petsc_options)
	~PETScLinearSolver ()
bool	solve (PETScMatrix &A, PETScVector &b, PETScVector &x)
PetscInt	getNumberOfIterations () const
	Get number of iterations.
double	getElapsedTime () const
	Get elapsed wall clock time.
bool	canSolveRectangular () const
	Get, if the solver can handle rectangular equation systems.
bool	willCompute (MathLib::LinearSolverBehaviour const) const
	Provided for compatibility with Eigen linear solvers.

Private Member Functions
bool	canSolverHandleRectangular (std::string_view ksp_type)

Private Attributes
KSP	solver_
	Solver type.
PC	pc_
	Preconditioner type.
bool	can_solve_rectangular_ = false
double	elapsed_ctime_ = 0.0
	Clock time.

Constructor & Destructor Documentation

◆ PETScLinearSolver()

MathLib::PETScLinearSolver::PETScLinearSolver	(	std::string const &	prefix,
		std::string const &	petsc_options )

Constructor

Parameters

prefix	Name used to distinguish the options in the command line for this solver. It can be the name of the PDE that owns an instance of this class.
petsc_options	PETSc options string which is passed to PETSc lib and inserted in the PETSc option database (see https://petsc.org/release/manualpages/Sys/PetscOptionsInsertString/).

Definition at line 10 of file PETScLinearSolver.cpp.

{
#if PETSC_VERSION_LT(3, 7, 0)
    PetscOptionsInsertString(petsc_options.c_str());
#else
    PetscOptionsInsertString(nullptr, petsc_options.c_str());
#endif
 
    PetscCallAbort(PETSC_COMM_WORLD, KSPCreate(PETSC_COMM_WORLD, &solver_));
 
    PetscCallAbort(PETSC_COMM_WORLD, KSPGetPC(solver_, &pc_));
 
    if (!prefix.empty())
    {
        PetscCallAbort(PETSC_COMM_WORLD,
                       KSPSetOptionsPrefix(solver_, prefix.c_str()));
    }
 
    PetscCallAbort(PETSC_COMM_WORLD,
                   KSPSetInitialGuessNonzero(solver_, PETSC_TRUE));
    PetscCallAbort(PETSC_COMM_WORLD,
                   KSPSetFromOptions(solver_));  // set run-time options
 
    KSPType ksp_type;
    PetscCallAbort(PETSC_COMM_WORLD, KSPGetType(solver_, &ksp_type));
    can_solve_rectangular_ = canSolverHandleRectangular(ksp_type);
}

References can_solve_rectangular_, canSolverHandleRectangular(), pc_, and solver_.

◆ ~PETScLinearSolver()

MathLib::PETScLinearSolver::~PETScLinearSolver ( )

inline

Definition at line 40 of file PETScLinearSolver.h.

    {
        PetscCallAbort(PETSC_COMM_WORLD, KSPDestroy(&solver_));
    }

References solver_.

Member Function Documentation

◆ canSolveRectangular()

bool MathLib::PETScLinearSolver::canSolveRectangular ( ) const

inline

Get, if the solver can handle rectangular equation systems.

Definition at line 59 of file PETScLinearSolver.h.

59{ return can_solve_rectangular_; }

References can_solve_rectangular_.

◆ canSolverHandleRectangular()

bool MathLib::PETScLinearSolver::canSolverHandleRectangular ( std::string_view ksp_type )

inlineprivate

Definition at line 69 of file PETScLinearSolver.h.

    {
        // List of KSP types that can handle rectangular matrices
        constexpr auto rectangular_solvers = std::to_array<std::string_view>(
            {KSPGMRES, KSPFGMRES, KSPBCGS, KSPBCGSL, KSPTFQMR, KSPCGS});
 
        return std::ranges::any_of(rectangular_solvers,
                                   [&](const auto& solver)
                                   { return solver == ksp_type; });
    }

Referenced by PETScLinearSolver().

◆ getElapsedTime()

double MathLib::PETScLinearSolver::getElapsedTime ( ) const

inline

Get elapsed wall clock time.

Definition at line 56 of file PETScLinearSolver.h.

56{ return elapsed_ctime_; }

MathLib::PETScLinearSolver::elapsed_ctime_

double elapsed_ctime_

Clock time.

Definition PETScLinearSolver.h:85

References elapsed_ctime_.

◆ getNumberOfIterations()

PetscInt MathLib::PETScLinearSolver::getNumberOfIterations ( ) const

inline

Get number of iterations.

Definition at line 48 of file PETScLinearSolver.h.

    {
        PetscInt its = 0;
        PetscCallAbort(PETSC_COMM_WORLD, KSPGetIterationNumber(solver_, &its));
        return its;
    }

References solver_.

◆ solve()

bool MathLib::PETScLinearSolver::solve	(	PETScMatrix &	A,
		PETScVector &	b,
		PETScVector &	x )

Definition at line 39 of file PETScLinearSolver.cpp.

{
    BaseLib::RunTime wtimer;
    wtimer.start();
 
// define TEST_MEM_PETSC
#ifdef TEST_MEM_PETSC
    PetscLogDouble mem1, mem2;
    PetscMemoryGetCurrentUsage(&mem1);
#endif
 
    KSPNormType norm_type;
    PetscCallAbort(PETSC_COMM_WORLD, KSPGetNormType(solver_, &norm_type));
    const char* ksp_type;
    const char* pc_type;
    PetscCallAbort(PETSC_COMM_WORLD, KSPGetType(solver_, &ksp_type));
    PetscCallAbort(PETSC_COMM_WORLD, PCGetType(pc_, &pc_type));
 
    PetscPrintf(PETSC_COMM_WORLD,
                "\n================================================");
    PetscPrintf(PETSC_COMM_WORLD,
                "\nLinear solver %s with %s preconditioner using %s", ksp_type,
                pc_type, KSPNormTypes[norm_type]);
 
    PetscCallAbort(PETSC_COMM_WORLD, KSPSetOperators(solver_, A.getRawMatrix(),
                                                     A.getRawMatrix()));
 
    PetscCallAbort(PETSC_COMM_WORLD,
                   KSPSolve(solver_, b.getRawVector(), x.getRawVector()));
 
    KSPConvergedReason reason;
    PetscCallAbort(PETSC_COMM_WORLD, KSPGetConvergedReason(solver_, &reason));
 
    bool converged = true;
    if (reason > 0)
    {
        PetscInt its;
        PetscCallAbort(PETSC_COMM_WORLD, KSPGetIterationNumber(solver_, &its));
        PetscPrintf(PETSC_COMM_WORLD, "\nconverged in %d iterations", its);
        switch (reason)
        {
            case KSP_CONVERGED_RTOL:
                PetscPrintf(PETSC_COMM_WORLD,
                            " (relative convergence criterion fulfilled).");
                break;
            case KSP_CONVERGED_ATOL:
                PetscPrintf(PETSC_COMM_WORLD,
                            " (absolute convergence criterion fulfilled).");
                break;
            default:
                PetscPrintf(PETSC_COMM_WORLD, ".");
        }
 
        PetscPrintf(PETSC_COMM_WORLD,
                    "\n================================================\n");
    }
    else if (reason == KSP_DIVERGED_ITS)
    {
        PetscPrintf(PETSC_COMM_WORLD,
                    "\nWarning: maximum number of iterations reached.\n");
    }
    else
    {
        converged = false;
        if (reason == KSP_DIVERGED_INDEFINITE_PC)
        {
            PetscPrintf(PETSC_COMM_WORLD,
                        "\nDivergence because of indefinite preconditioner,");
            PetscPrintf(PETSC_COMM_WORLD,
                        "\nTry to run again with "
                        "-pc_factor_shift_positive_definite option.\n");
        }
        else if (reason == KSP_DIVERGED_BREAKDOWN_BICG)
        {
            PetscPrintf(PETSC_COMM_WORLD,
                        "\nKSPBICG method was detected so the method could not "
                        "continue to enlarge the Krylov space.");
            PetscPrintf(PETSC_COMM_WORLD,
                        "\nTry to run again with another solver.\n");
        }
        else if (reason == KSP_DIVERGED_NONSYMMETRIC)
        {
            PetscPrintf(PETSC_COMM_WORLD,
                        "\nMatrix or preconditioner is unsymmetric but KSP "
                        "requires symmetric.\n");
        }
        else
        {
            PetscPrintf(PETSC_COMM_WORLD,
                        "\nDivergence detected, use command option "
                        "-ksp_monitor or -log_summary to check the details.\n");
        }
    }
 
#ifdef TEST_MEM_PETSC
    PetscMemoryGetCurrentUsage(&mem2);
    PetscPrintf(
        PETSC_COMM_WORLD,
        "###Memory usage by solver. Before: %f After: %f Increase: %d\n", mem1,
        mem2, (int)(mem2 - mem1));
#endif
 
    elapsed_ctime_ += wtimer.elapsed();
 
    return converged;
}

References BaseLib::RunTime::elapsed(), elapsed_ctime_, MathLib::PETScMatrix::getRawMatrix(), MathLib::PETScVector::getRawVector(), pc_, solver_, and BaseLib::RunTime::start().

◆ willCompute()

bool MathLib::PETScLinearSolver::willCompute ( MathLib::LinearSolverBehaviour const ) const

inline

Provided for compatibility with Eigen linear solvers.

Definition at line 62 of file PETScLinearSolver.h.

    {
        return true;
    }

Member Data Documentation

◆ can_solve_rectangular_

bool MathLib::PETScLinearSolver::can_solve_rectangular_ = false

private

Definition at line 83 of file PETScLinearSolver.h.

Referenced by PETScLinearSolver(), and canSolveRectangular().

◆ elapsed_ctime_

double MathLib::PETScLinearSolver::elapsed_ctime_ = 0.0

private

Clock time.

Definition at line 85 of file PETScLinearSolver.h.

Referenced by getElapsedTime(), and solve().

◆ pc_

PC MathLib::PETScLinearSolver::pc_

private

Preconditioner type.

Definition at line 81 of file PETScLinearSolver.h.

Referenced by PETScLinearSolver(), and solve().

◆ solver_

KSP MathLib::PETScLinearSolver::solver_

private

Solver type.

Definition at line 80 of file PETScLinearSolver.h.

Referenced by PETScLinearSolver(), ~PETScLinearSolver(), getNumberOfIterations(), and solve().

The documentation for this class was generated from the following files:

Detailed Description

Public Member Functions

Private Member Functions

Private Attributes

Constructor & Destructor Documentation

◆ PETScLinearSolver()

◆ ~PETScLinearSolver()

Member Function Documentation

◆ canSolveRectangular()

◆ canSolverHandleRectangular()

◆ getElapsedTime()

◆ getNumberOfIterations()

◆ solve()

◆ willCompute()

Member Data Documentation

◆ can_solve_rectangular_

◆ elapsed_ctime_

◆ pc_

◆ solver_