Loading...
Searching...
No Matches
AmgxPreconditioner.hpp
/*
Copyright 2024 SINTEF AS
Copyright 2024 Equinor ASA
This file is part of the Open Porous Media project (OPM).
OPM is free software: you can redistribute it and/or modify
it under the terms of the GNU General Public License as published by
the Free Software Foundation, either version 3 of the License, or
(at your option) any later version.
OPM is distributed in the hope that it will be useful,
but WITHOUT ANY WARRANTY; without even the implied warranty of
MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
GNU General Public License for more details.
You should have received a copy of the GNU General Public License
along with OPM. If not, see <http://www.gnu.org/licenses/>.
*/
#ifndef OPM_AMGX_PRECONDITIONER_HEADER_INCLUDED
#define OPM_AMGX_PRECONDITIONER_HEADER_INCLUDED
#include <opm/common/ErrorMacros.hpp>
#include <opm/common/TimingMacros.hpp>
#include <opm/simulators/linalg/PreconditionerWithUpdate.hpp>
#include <opm/simulators/linalg/PropertyTree.hpp>
#include <opm/simulators/linalg/gpuistl/AmgxInterface.hpp>
#include <dune/common/fmatrix.hh>
#include <dune/istl/bcrsmatrix.hh>
#include <amgx_c.h>
namespace Amgx
{
using AmgxInterface = Opm::gpuistl::AmgxInterface;
struct AmgxConfig {
int determinism_flag = 0;
int print_grid_stats = 0;
int print_solve_stats = 0;
std::string solver = "AMG";
std::string algorithm = "CLASSICAL";
std::string interpolator = "D2";
std::string selector = "PMIS";
std::string smoother = "BLOCK_JACOBI";
int presweeps = 3;
int postsweeps = 3;
double strength_threshold = 0.5;
int max_iters = 1;
explicit AmgxConfig(const Opm::PropertyTree& prm)
{
determinism_flag = prm.get<int>("determinism_flag", determinism_flag);
print_grid_stats = prm.get<int>("print_grid_stats", print_grid_stats);
print_solve_stats = prm.get<int>("print_solve_stats", print_solve_stats);
solver = prm.get<std::string>("solver", solver);
algorithm = prm.get<std::string>("algorithm", algorithm);
interpolator = prm.get<std::string>("interpolator", interpolator);
selector = prm.get<std::string>("selector", selector);
smoother = prm.get<std::string>("smoother", smoother);
presweeps = prm.get<int>("presweeps", presweeps);
postsweeps = prm.get<int>("postsweeps", postsweeps);
strength_threshold = prm.get<double>("strength_threshold", strength_threshold);
max_iters = prm.get<int>("max_iters", max_iters);
}
std::string toString() const {
return "config_version=2, "
"determinism_flag=" + std::to_string(determinism_flag) + ", "
"print_grid_stats=" + std::to_string(print_grid_stats) + ", "
"print_solve_stats=" + std::to_string(print_solve_stats) + ", "
"solver=" + solver + ", "
"algorithm=" + algorithm + ", "
"interpolator=" + interpolator + ", "
"selector=" + selector + ", "
"smoother=" + smoother + ", "
"presweeps=" + std::to_string(presweeps) + ", "
"postsweeps=" + std::to_string(postsweeps) + ", "
"strength_threshold=" + std::to_string(strength_threshold) + ", "
"max_iters=" + std::to_string(max_iters);
}
};
template <class M, class X, class Y>
class AmgxPreconditioner : public Dune::PreconditionerWithUpdate<X, Y>
{
public:
using matrix_type = M;
using matrix_field_type = typename M::field_type;
using domain_type = X;
using range_type = Y;
using vector_field_type = typename X::field_type;
static constexpr int block_size = 1;
AmgxPreconditioner(const M& A, const Opm::PropertyTree prm)
: A_(A)
, setup_frequency_(prm.get<int>("setup_frequency", 30))
, update_counter_(0)
{
OPM_TIMEBLOCK(prec_construct);
// Create configuration
AmgxConfig config(prm);
cfg_ = AmgxInterface::createConfig(config.toString());
rsrc_ = AmgxInterface::createResources(cfg_);
// Determine appropriate AMGX mode based on matrix and vector types
amgx_mode_ = AmgxInterface::determineAmgxMode<matrix_field_type, vector_field_type>();
// Create solver, matrix, and vector handles, given the mode
solver_ = AmgxInterface::createSolver(rsrc_, amgx_mode_, cfg_);
A_amgx_ = AmgxInterface::createMatrix(rsrc_, amgx_mode_);
x_amgx_ = AmgxInterface::createVector(rsrc_, amgx_mode_);
b_amgx_ = AmgxInterface::createVector(rsrc_, amgx_mode_);
// Initialize matrix structure and values
AmgxInterface::initializeMatrix(A_, A_amgx_);
// Initialize vectors with proper dimensions
const int N = A_.N();
AmgxInterface::initializeVector(N, block_size, x_amgx_);
AmgxInterface::initializeVector(N, block_size, b_amgx_);
// Perform initial update
update();
}
~AmgxPreconditioner()
{
AmgxInterface::destroySolver(solver_);
AmgxInterface::destroyVector(x_amgx_);
AmgxInterface::destroyVector(b_amgx_);
AmgxInterface::destroyMatrix(A_amgx_, A_);
AmgxInterface::destroyResources(rsrc_);
AmgxInterface::destroyConfig(cfg_);
}
void pre(X& /*v*/, Y& /*d*/) override
{
}
void apply(X& v, const Y& d) override
{
OPM_TIMEBLOCK(prec_apply);
// Transfer vectors to AMGX
AmgxInterface::transferVectorToAmgx(v, x_amgx_);
AmgxInterface::transferVectorToAmgx(d, b_amgx_);
// Apply preconditioner
AMGX_SAFE_CALL(AMGX_solver_solve(solver_, b_amgx_, x_amgx_));
// Transfer result back
AmgxInterface::transferVectorFromAmgx(x_amgx_, v);
}
void post(X& /*v*/) override
{
}
void update() override
{
OPM_TIMEBLOCK(prec_update);
AmgxInterface::updateMatrixValues(A_, A_amgx_);
// Setup or resetup the solver based on update counter
if (update_counter_ == 0) {
AMGX_SAFE_CALL(AMGX_solver_setup(solver_, A_amgx_));
} else {
AMGX_SAFE_CALL(AMGX_solver_resetup(solver_, A_amgx_));
}
++update_counter_;
if (update_counter_ >= setup_frequency_) {
update_counter_ = 0;
}
}
Dune::SolverCategory::Category category() const override
{
return Dune::SolverCategory::sequential;
}
bool hasPerfectUpdate() const override
{
// The AMG hierarchy of the Amgx preconditioner can depend on the values of the matrix, so it must be recreated
// when the matrix values change, at a given frequency. Since this is handled internally, we return true.
return true;
}
private:
// Reference to the input matrix
const M& A_;
// AMGX handles
AMGX_config_handle cfg_ = nullptr;
AMGX_resources_handle rsrc_ = nullptr;
AMGX_solver_handle solver_ = nullptr;
AMGX_matrix_handle A_amgx_ = nullptr;
AMGX_vector_handle x_amgx_ = nullptr;
AMGX_vector_handle b_amgx_ = nullptr;
AMGX_Mode amgx_mode_;
// Frequency of updating the AMG hierarchy
int setup_frequency_;
// Counter for setup updates
int update_counter_;
};
} // namespace Amgx
#endif // OPM_AMGX_PRECONDITIONER_HEADER_INCLUDED
Amgx::AmgxPreconditioner::hasPerfectUpdate
bool hasPerfectUpdate() const override
Checks if the preconditioner has a perfect update.
Definition AmgxPreconditioner.hpp:265
Amgx::AmgxPreconditioner::range_type
Y range_type
The range type of the preconditioner.
Definition AmgxPreconditioner.hpp:116
Amgx::AmgxPreconditioner::update
void update() override
Updates the preconditioner with the current matrix values.
Definition AmgxPreconditioner.hpp:231
Amgx::AmgxPreconditioner::apply
void apply(X &v, const Y &d) override
Applies the preconditioner to a vector.
Definition AmgxPreconditioner.hpp:200
Amgx::AmgxPreconditioner::matrix_field_type
typename M::field_type matrix_field_type
The field type of the matrix.
Definition AmgxPreconditioner.hpp:112
Amgx::AmgxPreconditioner::AmgxPreconditioner
AmgxPreconditioner(const M &A, const Opm::PropertyTree prm)
Constructor for AmgxPreconditioner.
Definition AmgxPreconditioner.hpp:130
Amgx::AmgxPreconditioner::category
Dune::SolverCategory::Category category() const override
Returns the solver category.
Definition AmgxPreconditioner.hpp:255
Amgx::AmgxPreconditioner::pre
void pre(X &, Y &) override
Pre-processing step before applying the preconditioner.
Definition AmgxPreconditioner.hpp:186
Amgx::AmgxPreconditioner::domain_type
X domain_type
The domain type of the preconditioner.
Definition AmgxPreconditioner.hpp:114
Amgx::AmgxPreconditioner::matrix_type
M matrix_type
The matrix type the preconditioner is for.
Definition AmgxPreconditioner.hpp:110
Amgx::AmgxPreconditioner::vector_field_type
typename X::field_type vector_field_type
The field type of the vectors.
Definition AmgxPreconditioner.hpp:118
Amgx::AmgxPreconditioner::~AmgxPreconditioner
~AmgxPreconditioner()
Destructor for AmgxPreconditioner.
Definition AmgxPreconditioner.hpp:168
Amgx::AmgxPreconditioner::post
void post(X &) override
Post-processing step after applying the preconditioner.
Definition AmgxPreconditioner.hpp:222
Dune::PreconditionerWithUpdate
Interface class adding the update() method to the preconditioner interface.
Definition PreconditionerWithUpdate.hpp:32
Opm::PropertyTree
Hierarchical collection of key/value pairs.
Definition PropertyTree.hpp:39
Opm::PropertyTree::get
T get(const std::string &key) const
Retrieve property value given hierarchical property key.
Definition PropertyTree.cpp:59
Opm::gpuistl::AmgxInterface::determineAmgxMode
static AMGX_Mode determineAmgxMode()
Determine the appropriate AMGX mode based on matrix and vector field types.
Definition AmgxInterface.hpp:613
Opm::gpuistl::AmgxInterface::destroyResources
static void destroyResources(AMGX_resources_handle resources)
Destroy an AMGX resources handle.
Definition AmgxInterface.hpp:244
Opm::gpuistl::AmgxInterface::createResources
static AMGX_resources_handle createResources(AMGX_config_handle config)
Create AMGX resources from a config.
Definition AmgxInterface.hpp:167
Opm::gpuistl::AmgxInterface::destroyVector
static void destroyVector(AMGX_vector_handle vector)
Destroy an AMGX vector handle.
Definition AmgxInterface.hpp:292
Opm::gpuistl::AmgxInterface::createSolver
static AMGX_solver_handle createSolver(AMGX_resources_handle resources, AMGX_Mode mode, AMGX_config_handle config)
Create an AMGX solver.
Definition AmgxInterface.hpp:183
Opm::gpuistl::AmgxInterface::updateMatrixValues
static void updateMatrixValues(const MatrixType &matrix, AMGX_matrix_handle amgx_matrix)
Update matrix values in AMGX.
Definition AmgxInterface.hpp:590
Opm::gpuistl::AmgxInterface::initializeVector
static void initializeVector(int N, int block_size, AMGX_vector_handle amgx_vector)
Initialize an AMGX vector with zeros.
Definition AmgxInterface.hpp:574
Opm::gpuistl::AmgxInterface::destroyConfig
static void destroyConfig(AMGX_config_handle config)
Destroy an AMGX config handle.
Definition AmgxInterface.hpp:226
Opm::gpuistl::AmgxInterface::createConfig
static AMGX_config_handle createConfig(const std::string &config_string)
Create an AMGX config handle from a configuration string.
Definition AmgxInterface.hpp:153
Opm::gpuistl::AmgxInterface::destroySolver
static void destroySolver(AMGX_solver_handle solver)
Destroy an AMGX solver handle.
Definition AmgxInterface.hpp:257
Opm::gpuistl::AmgxInterface::createMatrix
static AMGX_matrix_handle createMatrix(AMGX_resources_handle resources, AMGX_Mode mode)
Create an AMGX matrix.
Definition AmgxInterface.hpp:198
Opm::gpuistl::AmgxInterface::initializeMatrix
static void initializeMatrix(const MatrixType &matrix, AMGX_matrix_handle amgx_matrix)
Initialize an AMGX matrix from any matrix type (CPU or GPU).
Definition AmgxInterface.hpp:515
Opm::gpuistl::AmgxInterface::createVector
static AMGX_vector_handle createVector(AMGX_resources_handle resources, AMGX_Mode mode)
Create an AMGX vector.
Definition AmgxInterface.hpp:213
Opm::gpuistl::AmgxInterface::destroyMatrix
static void destroyMatrix(AMGX_matrix_handle amgx_matrix, const MatrixType &matrix)
Destroy an AMGX matrix handle.
Definition AmgxInterface.hpp:272
Opm::gpuistl::AmgxInterface::transferVectorToAmgx
static void transferVectorToAmgx(const VectorType &vec, AMGX_vector_handle amgx_vec)
Transfer vector to AMGX from any vector type (CPU or GPU).
Definition AmgxInterface.hpp:370
Opm::gpuistl::AmgxInterface::transferVectorFromAmgx
static void transferVectorFromAmgx(AMGX_vector_handle amgx_vec, VectorType &vec)
Transfer vector from AMGX to any vector type (CPU or GPU).
Definition AmgxInterface.hpp:396
Amgx::AmgxConfig
Configuration structure for AMGX parameters.
Definition AmgxPreconditioner.hpp:45