OilfieldToolsNet/opm-simulators
4
0 Fork 0
mirror of https://github.com/OPM/opm-simulators.git synced 2026-10-25 15:10:09 -05:00
Code Issues Packages Projects Releases Wiki Activity

Removing the ISTLSolverGPUISTLImplementation

Browse Source
This commit is contained in:
Kjetil Olsen Lye committed 2025-06-16 11:41:32 +02:00
1 parent f30036a96c
commit 12bd937a26
4 files changed
+69 -271
Download Patch File Download Diff File Expand all files Collapse all files

No files matched your search

CMakeLists_files.cmake
-3
View File
@@ -306,8 +306,6 @@ if (HAVE_CUDA)
ADD_CUDA_OR_HIP_FILE(MAIN_SOURCE_FILES opm/simulators/linalg detail/FlexibleSolverWrapper.cpp)
ADD_CUDA_OR_HIP_FILE(MAIN_SOURCE_FILES opm/simulators/linalg FlexibleSolver_gpu_instantiate.cpp)
ADD_CUDA_OR_HIP_FILE(MAIN_SOURCE_FILES opm/simulators/linalg PreconditionerFactory_gpu_instantiate.cpp)
ADD_CUDA_OR_HIP_FILE(MAIN_SOURCE_FILES opm/simulators/linalg detail/ISTLSolverGPUISTLImplementation.cpp)
# HEADERS
@@ -356,7 +354,6 @@ if (HAVE_CUDA)
ADD_CUDA_OR_HIP_FILE(PUBLIC_HEADER_FILES opm/simulators/linalg PreconditionerCPUMatrixToGPUMatrix.hpp)
ADD_CUDA_OR_HIP_FILE(PUBLIC_HEADER_FILES opm/simulators/linalg ISTLSolverGPUISTL.hpp)
ADD_CUDA_OR_HIP_FILE(PUBLIC_HEADER_FILES opm/simulators/linalg detail/FlexibleSolverWrapper.hpp)
ADD_CUDA_OR_HIP_FILE(PUBLIC_HEADER_FILES opm/simulators/linalg detail/ISTLSolverGPUISTLImplementation.hpp)
if(MPI_FOUND)
opm/simulators/linalg/gpuistl/ISTLSolverGPUISTL.hpp
+69 -12
View File
@@ -22,8 +22,6 @@
#include <opm/simulators/linalg/AbstractISTLSolver.hpp>
#include <opm/simulators/linalg/ISTLSolver.hpp>
#include <opm/simulators/linalg/gpuistl/detail/ISTLSolverGPUISTLImplementation.hpp>
#if USE_HIP
#include <opm/simulators/linalg/gpuistl_hip/GpuVector.hpp>
#include <opm/simulators/linalg/gpuistl_hip/GpuSparseMatrix.hpp>
@@ -32,6 +30,8 @@
#include <opm/simulators/linalg/gpuistl/GpuSparseMatrix.hpp>
#endif
#include <opm/simulators/linalg/gpuistl/detail/FlexibleSolverWrapper.hpp>
namespace Opm::gpuistl
{
@@ -44,6 +44,12 @@ public:
using Vector = GetPropType<TypeTag, Properties::GlobalEqVector>;
using Simulator = GetPropType<TypeTag, Properties::Simulator>;
using Matrix = typename SparseMatrixAdapter::IstlMatrix;
using real_type = typename Vector::field_type;
using GPUMatrix = Opm::gpuistl::GpuSparseMatrix<real_type>;
using GPUVector = Opm::gpuistl::GpuVector<real_type>;
constexpr static std::size_t pressureIndex = GetPropType<TypeTag, Properties::Indices>::pressureSwitchIdx;
@@ -53,8 +59,7 @@ public:
using CommunicationType = Dune::Communication<int>;
#endif
using SolverType = Opm::gpuistl::detail::FlexibleSolverWrapper<GPUMatrix, GPUVector, CommunicationType>;
/// Construct a system solver.
/// \param[in] simulator The opm-models simulator object
@@ -66,14 +71,12 @@ public:
ISTLSolverGPUISTL(const Simulator& simulator,
const FlowLinearSolverParameters& parameters,
bool forceSerial = false)
: m_parameters(parameters)
: m_parameters(parameters), m_forceSerial(forceSerial)
{
m_parameters.init(simulator.vanguard().eclState().getSimulationConfig().useCPR());
m_propertyTree = setupPropertyTree(m_parameters,
Parameters::IsSet<Parameters::LinearSolverMaxIter>(),
Parameters::IsSet<Parameters::LinearSolverReduction>());
m_implementation = std::make_unique<detail::ISTLSolverGPUISTLImplementation<Matrix, Vector, CommunicationType>>(
m_parameters, m_propertyTree, forceSerial, pressureIndex);
}
/// Construct a system solver.
@@ -110,7 +113,8 @@ public:
void prepare(const Matrix& M, Vector& b) override
{
try {
m_implementation->prepare(M, b);
updateMatrix(M);
updateRhs(b);
}
OPM_CATCH_AND_RETHROW_AS_CRITICAL_ERROR("This is likely due to a faulty linear solver JSON specification. "
"Check for errors related to missing nodes.");
@@ -123,7 +127,10 @@ public:
void getResidual(Vector& b) const override
{
m_implementation->getResidual(b);
if (!m_rhs) {
OPM_THROW(std::runtime_error, "m_rhs not initialized, prepare(matrix, rhs); needs to be called");
}
m_rhs->copyToHost(b);
}
void setMatrix(const SparseMatrixAdapter&) override
@@ -135,8 +142,28 @@ public:
{
// TODO: Write matrix to disk if needed
Dune::InverseOperatorResult result;
m_implementation->apply(x, result);
if (!m_matrix) {
OPM_THROW(std::runtime_error, "m_matrix not initialized, prepare(matrix, rhs); needs to be called");
}
if (!m_rhs) {
OPM_THROW(std::runtime_error, "m_rhs not initialized, prepare(matrix, rhs); needs to be called");
}
if (!m_gpuSolver) {
OPM_THROW(std::runtime_error, "m_gpuFlexibleSolver not initialized, prepare(matrix, rhs); needs to be called");
}
if (!m_x) {
m_x = std::make_unique<GpuVector<real_type>>(x);
} else {
m_x->copyFromHost(x);
}
m_gpuSolver->apply(*m_x, *m_rhs, result);
m_x->copyToHost(x);
++m_solveCount;
m_lastSeenIterations = result.iterations;
if (!result.converged) {
if (result.reduction < m_parameters.relaxed_linear_solver_reduction_) {
@@ -153,7 +180,6 @@ public:
OPM_THROW_NOLOG(NumericalProblem, msg);
}
return result.converged;
}
@@ -174,12 +200,43 @@ public:
}
private:
void updateMatrix(const Matrix& M) {
if (!m_matrix) {
m_matrix.reset(new auto(GPUMatrix::fromMatrix(M)));
std::function<GPUVector()> weightsCalculator = {};
const bool parallel = false;
m_gpuSolver = std::make_unique<SolverType>(
*m_matrix, parallel, m_propertyTree, pressureIndex, weightsCalculator, m_forceSerial, nullptr);
} else {
m_matrix->updateNonzeroValues(M);
}
m_gpuSolver->update();
}
void updateRhs(const Vector& b) {
if (!m_rhs) {
m_rhs = std::make_unique<GPUVector>(b);
} else {
m_rhs->copyFromHost(b);
}
}
FlowLinearSolverParameters m_parameters;
PropertyTree m_propertyTree;
std::unique_ptr<detail::ISTLSolverGPUISTLImplementation<Matrix, Vector, CommunicationType>> m_implementation;
int m_lastSeenIterations = 0;
int m_solveCount = 0;
std::unique_ptr<GPUMatrix> m_matrix;
std::unique_ptr<SolverType> m_gpuSolver;
std::unique_ptr<GPUVector> m_rhs;
std::unique_ptr<GPUVector> m_x;
const bool m_forceSerial;
};
} // namespace Opm::gpuistl
opm/simulators/linalg/gpuistl/detail/ISTLSolverGPUISTLImplementation.cpp
-156
View File
@@ -1,156 +0,0 @@
/*
Copyright 2025 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/>.
*/
#include "config.h"
#include "dune/common/fvector.hh"
#include "dune/istl/bvector.hh"
#include <dune/common/fmatrix.hh>
#include <dune/istl/bcrsmatrix.hh>
#include <opm/simulators/linalg/gpuistl/GpuSparseMatrix.hpp>
#include <opm/simulators/linalg/gpuistl/GpuVector.hpp>
#include <opm/simulators/linalg/gpuistl/detail/ISTLSolverGPUISTLImplementation.hpp>
#include <opm/simulators/linalg/matrixblock.hh>
namespace Opm::gpuistl::detail
{
template <class CpuMatrix, class CpuVector, class Comm>
ISTLSolverGPUISTLImplementation<CpuMatrix, CpuVector, Comm>::ISTLSolverGPUISTLImplementation(
const FlowLinearSolverParameters& parameters,
const PropertyTree& propertyTree,
bool forceSerial,
std::size_t pressureIndex)
: m_parameters(parameters)
, m_propertyTree(propertyTree)
, m_forceSerial(forceSerial)
, m_pressureIndex(pressureIndex)
{
}
template <class CpuMatrix, class CpuVector, class Comm>
void
ISTLSolverGPUISTLImplementation<CpuMatrix, CpuVector, Comm>::prepare(const CpuMatrix& M, CpuVector& b)
{
updateMatrix(M);
updateRhs(b);
}
template <class CpuMatrix, class CpuVector, class Comm>
void
ISTLSolverGPUISTLImplementation<CpuMatrix, CpuVector, Comm>::apply(CpuVector& x, Dune::InverseOperatorResult& result)
{
if (!m_matrix) {
OPM_THROW(std::runtime_error, "m_matrix not initialized, prepare(matrix, rhs); needs to be called");
}
if (!m_rhs) {
OPM_THROW(std::runtime_error, "m_rhs not initialized, prepare(matrix, rhs); needs to be called");
}
if (!m_gpuSolver) {
OPM_THROW(std::runtime_error, "m_gpuFlexibleSolver not initialized, prepare(matrix, rhs); needs to be called");
}
if (!m_x) {
m_x = std::make_unique<GpuVector<real_type>>(x);
} else {
m_x->copyFromHost(x);
}
m_gpuSolver->apply(*m_x, *m_rhs, result);
m_x->copyToHost(x);
// The rest is handled "upstream"
}
template <class CpuMatrix, class CpuVector, class Comm>
void
ISTLSolverGPUISTLImplementation<CpuMatrix, CpuVector, Comm>::getResidual(CpuVector& b) const
{
if (!m_rhs) {
OPM_THROW(std::runtime_error, "m_rhs not initialized, prepare(matrix, rhs); needs to be called");
}
m_rhs->copyToHost(b);
}
template <class CpuMatrix, class CpuVector, class Comm>
void
ISTLSolverGPUISTLImplementation<CpuMatrix, CpuVector, Comm>::updateMatrix(const CpuMatrix& M)
{
if (!m_matrix) {
m_matrix.reset(new auto(GPUMatrix::fromMatrix(M)));
std::function<GPUVector()> weightsCalculator = {};
const bool parallel = false;
m_gpuSolver = std::make_unique<SolverType>(
*m_matrix, parallel, m_propertyTree, m_pressureIndex, weightsCalculator, m_forceSerial, nullptr);
} else {
m_matrix->updateNonzeroValues(M);
}
m_gpuSolver->update();
}
template <class CpuMatrix, class CpuVector, class Comm>
void
ISTLSolverGPUISTLImplementation<CpuMatrix, CpuVector, Comm>::updateRhs(const CpuVector& b)
{
if (!m_rhs) {
m_rhs = std::make_unique<GPUVector>(b);
} else {
m_rhs->copyFromHost(b);
}
}
template <class CpuMatrix, class CpuVector, class Comm>
ISTLSolverGPUISTLImplementation<CpuMatrix, CpuVector, Comm>::~ISTLSolverGPUISTLImplementation() = default;
} // namespace Opm::gpuistl::detail
#if HAVE_MPI
// Parallel communicator and operators.
using Comm = Dune::OwnerOverlapCopyCommunication<int, int>;
#else
using Comm = Dune::Communication<int>;
#endif
#define INSTANTIATE_ISTLSOLVER_GPUISTL_DUNE(realtype, blockdim) \
template class ::Opm::gpuistl::detail::ISTLSolverGPUISTLImplementation< \
::Dune::BCRSMatrix<Dune::FieldMatrix<realtype, blockdim, blockdim>>, \
::Dune::BlockVector<::Dune::FieldVector<realtype, blockdim>>, \
Comm>; \
template class ::Opm::gpuistl::detail::ISTLSolverGPUISTLImplementation< \
::Dune::BCRSMatrix<Opm::MatrixBlock<realtype, blockdim, blockdim>>, \
::Dune::BlockVector<::Dune::FieldVector<realtype, blockdim>>, \
Comm>
INSTANTIATE_ISTLSOLVER_GPUISTL_DUNE(double, 1);
INSTANTIATE_ISTLSOLVER_GPUISTL_DUNE(double, 2);
INSTANTIATE_ISTLSOLVER_GPUISTL_DUNE(double, 3);
INSTANTIATE_ISTLSOLVER_GPUISTL_DUNE(double, 4);
INSTANTIATE_ISTLSOLVER_GPUISTL_DUNE(double, 5);
INSTANTIATE_ISTLSOLVER_GPUISTL_DUNE(double, 6);
#if FLOW_INSTANTIATE_FLOAT
INSTANTIATE_ISTLSOLVER_GPUISTL_DUNE(float, 1);
INSTANTIATE_ISTLSOLVER_GPUISTL_DUNE(float, 2);
INSTANTIATE_ISTLSOLVER_GPUISTL_DUNE(float, 3);
INSTANTIATE_ISTLSOLVER_GPUISTL_DUNE(float, 4);
INSTANTIATE_ISTLSOLVER_GPUISTL_DUNE(float, 5);
INSTANTIATE_ISTLSOLVER_GPUISTL_DUNE(float, 6);
#endif
opm/simulators/linalg/gpuistl/detail/ISTLSolverGPUISTLImplementation.hpp
-100
View File
@@ -1,100 +0,0 @@
/*
Copyright 2025 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_ISTLSOLVERGPUISTLIMPLEMENTATION_HEADER_INCLUDED
#define OPM_ISTLSOLVERGPUISTLIMPLEMENTATION_HEADER_INCLUDED
#include "dune/istl/solver.hh"
#include <memory>
#include <dune/istl/operators.hh>
#include <opm/simulators/linalg/FlowLinearSolverParameters.hpp>
#include <opm/simulators/linalg/PropertyTree.hpp>
#include <opm/simulators/linalg/gpuistl/detail/FlexibleSolverWrapper.hpp>
namespace Opm::gpuistl
{
// We need to forward declare these classes, but not include the headers,
// to be able to include this header in main files.
//
// TODO: Fix hipification so that we can include the headers here.
// Once hipification is fixed, this class can be integrated into ISTLSolverGPUISTL.
template <class T>
class GpuSparseMatrix;
template <class T>
class GpuVector;
} // namespace Opm::gpuistl
namespace Opm::gpuistl::detail
{
/**
* This is an implementation detail class for the ISTLSolverGPUISTL class.
* It is not intended to be used directly.
*
* In short, this class just exists until we can make the hipification
* a bit more robust.
*/
template <class CpuMatrix, class CpuVector, class Comm>
class ISTLSolverGPUISTLImplementation
{
public:
using real_type = typename CpuVector::field_type;
using GPUMatrix = Opm::gpuistl::GpuSparseMatrix<real_type>;
using GPUVector = Opm::gpuistl::GpuVector<real_type>;
using SolverType = Opm::gpuistl::detail::FlexibleSolverWrapper<GPUMatrix, GPUVector, Comm>;
ISTLSolverGPUISTLImplementation(const FlowLinearSolverParameters& parameters,
const PropertyTree& propertyTree,
bool forceSerial,
std::size_t pressureIndex);
ISTLSolverGPUISTLImplementation() = delete;
ISTLSolverGPUISTLImplementation(const ISTLSolverGPUISTLImplementation&) = delete;
ISTLSolverGPUISTLImplementation& operator=(const ISTLSolverGPUISTLImplementation&) = delete;
~ISTLSolverGPUISTLImplementation();
void prepare(const CpuMatrix& M, CpuVector& b);
void apply(CpuVector& x, Dune::InverseOperatorResult& result);
void getResidual(CpuVector& b) const;
private:
void updateMatrix(const CpuMatrix& M);
void updateRhs(const CpuVector& b);
std::unique_ptr<GPUMatrix> m_matrix;
std::unique_ptr<SolverType> m_gpuSolver;
std::unique_ptr<GPUVector> m_rhs;
std::unique_ptr<GPUVector> m_x;
const FlowLinearSolverParameters& m_parameters;
const PropertyTree& m_propertyTree;
const bool m_forceSerial;
const std::size_t m_pressureIndex;
};
} // namespace Opm::gpuistl::detail
#endif // OPM_ISTLSOLVERGPUISTLIMPLEMENTATION_HEADER_INCLUDED