OilfieldToolsNet/opm-simulators
4
0
Fork 0
mirror of https://github.com/OPM/opm-simulators.git synced 2025-02-25 18:55:30 -06:00
Code Issues Packages Projects Releases Wiki Activity

changed: put parts of ISTLSolverEbos in separate compile unit

Browse Source 
in particular this means simulator objects are now not dependent
on FlexibleSolver.hpp and BdaBridge.hpp
This commit is contained in:
Arne Morten Kvarving 2022-08-17 10:37:00 +02:00
parent 3fecc92956
commit 8e2b385f4d
5 changed files with 566 additions and 267 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

1
CMakeLists_files.cmake
View File

@ -51,6 +51,7 @@ list (APPEND MAIN_SOURCE_FILES
opm/simulators/linalg/FlexibleSolver4.cpp
opm/simulators/linalg/FlexibleSolver5.cpp
opm/simulators/linalg/FlexibleSolver6.cpp
opm/simulators/linalg/ISTLSolverEbos.cpp
opm/simulators/linalg/MILU.cpp
opm/simulators/linalg/ParallelIstlInformation.cpp
opm/simulators/linalg/ParallelOverlappingILU0.cpp

1
ebos/eclgenerictracermodel.cc
View File

@ -33,7 +33,6 @@
#include <opm/input/eclipse/EclipseState/EclipseState.hpp>
#include <opm/input/eclipse/EclipseState/Runspec.hpp>
#include <opm/input/eclipse/EclipseState/Tables/TracerVdTable.hpp>
#include <opm/simulators/linalg/FlexibleSolver.hpp>
#include <dune/istl/operators.hh>
#include <dune/istl/solvers.hh>

371
opm/simulators/linalg/ISTLSolverEbos.cpp Normal file
View File

@ -0,0 +1,371 @@
/*
Copyright 2016 IRIS AS
Copyright 2019, 2020 Equinor ASA
Copyright 2020 SINTEF Digital, Mathematics and Cybernetics
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 <opm/simulators/linalg/ISTLSolverEbos.hpp>
#include <dune/istl/schwarz.hh>
#include <opm/grid/CpGrid.hpp>
#include <opm/simulators/linalg/FlexibleSolver.hpp>
#include <opm/simulators/linalg/ParallelIstlInformation.hpp>
#include <opm/simulators/utils/ParallelCommunication.hpp>
#if HAVE_CUDA || HAVE_OPENCL || HAVE_FPGA || HAVE_AMGCL
#include <opm/simulators/linalg/bda/BdaBridge.hpp>
#include <opm/simulators/linalg/bda/WellContributions.hpp>
#include <iostream>
#endif
namespace Opm {
namespace detail {
#ifdef HAVE_MPI
void copyParValues(std::any& parallelInformation, size_t size,
Dune::OwnerOverlapCopyCommunication<int,int>& comm)
{
if (parallelInformation.type() == typeid(ParallelISTLInformation)) {
const ParallelISTLInformation* parinfo = std::any_cast<ParallelISTLInformation>(&parallelInformation);
assert(parinfo);
parinfo->copyValuesTo(comm.indexSet(), comm.remoteIndices(), size, 1);
}
}
#endif
template<class Matrix>
void makeOverlapRowsInvalid(Matrix& matrix,
const std::vector<int>& overlapRows)
{
//value to set on diagonal
const int numEq = Matrix::block_type::rows;
typename Matrix::block_type diag_block(0.0);
for (int eq = 0; eq < numEq; ++eq)
diag_block[eq][eq] = 1.0;
//loop over precalculated overlap rows and columns
for (const auto row : overlapRows)
{
// Zero out row.
matrix[row] = 0.0;
//diagonal block set to diag(1.0).
matrix[row][row] = diag_block;
}
}
/// Return an appropriate weight function if a cpr preconditioner is asked for.
template<class Vector, class Matrix>
std::function<Vector()> getWeightsCalculator(const PropertyTree& prm,
const Matrix& matrix,
size_t pressureIndex,
std::function<Vector()> trueFunc)
{
std::function<Vector()> weightsCalculator;
using namespace std::string_literals;
auto preconditionerType = prm.get("preconditioner.type"s, "cpr"s);
if (preconditionerType == "cpr" || preconditionerType == "cprt"
|| preconditionerType == "cprw" || preconditionerType == "cprwt") {
const bool transpose = preconditionerType == "cprt" || preconditionerType == "cprwt";
const auto weightsType = prm.get("preconditioner.weight_type"s, "quasiimpes"s);
if (weightsType == "quasiimpes") {
// weights will be created as default in the solver
// assignment p = pressureIndex prevent compiler warning about
// capturing variable with non-automatic storage duration
weightsCalculator = [matrix, transpose, pressureIndex]() {
return Amg::getQuasiImpesWeights<Matrix, Vector>(matrix,
pressureIndex,
transpose);
};
} else if (weightsType == "trueimpes") {
weightsCalculator = trueFunc;
} else {
OPM_THROW(std::invalid_argument,
"Weights type " << weightsType << "not implemented for cpr."
<< " Please use quasiimpes or trueimpes.");
}
}
return weightsCalculator;
}
template<class Matrix, class Vector, class Comm>
void FlexibleSolverInfo<Matrix,Vector,Comm>::create(const Matrix& matrix,
bool parallel,
const PropertyTree& prm,
size_t pressureIndex,
std::function<Vector()> trueFunc,
[[maybe_unused]] Comm& comm)
{
std::function<Vector()> weightsCalculator =
getWeightsCalculator<Vector>(prm, matrix, pressureIndex, trueFunc);
if (parallel) {
#if HAVE_MPI
if (!wellOperator_) {
using ParOperatorType = Dune::OverlappingSchwarzOperator<Matrix, Vector, Vector, Comm>;
auto pop = std::make_unique<ParOperatorType>(matrix, comm);
using FlexibleSolverType = Dune::FlexibleSolver<ParOperatorType>;
auto sol = std::make_unique<FlexibleSolverType>(*pop, comm, prm,
weightsCalculator,
pressureIndex);
this->pre_ = &sol->preconditioner();
this->op_ = std::move(pop);
this->solver_ = std::move(sol);
} else {
using ParOperatorType = WellModelGhostLastMatrixAdapter<Matrix, Vector, Vector, true>;
auto pop = std::make_unique<ParOperatorType>(matrix, *wellOperator_,
interiorCellNum_);
using FlexibleSolverType = Dune::FlexibleSolver<ParOperatorType>;
auto sol = std::make_unique<FlexibleSolverType>(*pop, comm, prm,
weightsCalculator,
pressureIndex);
this->pre_ = &sol->preconditioner();
this->op_ = std::move(pop);
this->solver_ = std::move(sol);
}
#endif
} else {
if (!wellOperator_) {
using SeqOperatorType = Dune::MatrixAdapter<Matrix, Vector, Vector>;
auto sop = std::make_unique<SeqOperatorType>(matrix);
using FlexibleSolverType = Dune::FlexibleSolver<SeqOperatorType>;
auto sol = std::make_unique<FlexibleSolverType>(*sop, prm,
weightsCalculator,
pressureIndex);
this->pre_ = &sol->preconditioner();
this->op_ = std::move(sop);
this->solver_ = std::move(sol);
} else {
using SeqOperatorType = WellModelMatrixAdapter<Matrix, Vector, Vector, false>;
auto sop = std::make_unique<SeqOperatorType>(matrix, *wellOperator_);
using FlexibleSolverType = Dune::FlexibleSolver<SeqOperatorType>;
auto sol = std::make_unique<FlexibleSolverType>(*sop, prm,
weightsCalculator,
pressureIndex);
this->pre_ = &sol->preconditioner();
this->op_ = std::move(sop);
this->solver_ = std::move(sol);
}
}
}
#if HAVE_CUDA || HAVE_OPENCL || HAVE_FPGA || HAVE_AMGCL
template<class Matrix, class Vector>
BdaSolverInfo<Matrix,Vector>::
BdaSolverInfo(const std::string& accelerator_mode,
const std::string& fpga_bitstream,
const int linear_solver_verbosity,
const int maxit,
const double tolerance,
const int platformID,
const int deviceID,
const std::string& opencl_ilu_reorder,
const std::string& linsolver)
: bridge_(std::make_unique<Bridge>(accelerator_mode, fpga_bitstream,
linear_solver_verbosity, maxit,
tolerance, platformID, deviceID,
opencl_ilu_reorder, linsolver))
, accelerator_mode_(accelerator_mode)
{}
template<class Matrix, class Vector>
BdaSolverInfo<Matrix,Vector>::~BdaSolverInfo() = default;
template<class Matrix, class Vector>
template<class Grid>
void BdaSolverInfo<Matrix,Vector>::
prepare(const Grid& grid,
const Dune::CartesianIndexMapper<Grid>& cartMapper,
const std::vector<Well>& wellsForConn,
const std::vector<int>& cellPartition,
const size_t nonzeroes,
const bool useWellConn)
{
if (numJacobiBlocks_ > 1) {
detail::setWellConnections(grid, cartMapper, wellsForConn,
useWellConn,
wellConnectionsGraph_,
numJacobiBlocks_);
std::cout << "Create block-Jacobi pattern" << std::endl;
this->blockJacobiAdjacency(grid, cellPartition, nonzeroes);
}
}
template<class Matrix, class Vector>
bool BdaSolverInfo<Matrix,Vector>::
apply(Vector& rhs,
const bool useWellConn,
WellContribFunc getContribs,
const int rank,
Matrix& matrix,
Vector& x,
Dune::InverseOperatorResult& result)
{
bool use_gpu = bridge_->getUseGpu();
bool use_fpga = bridge_->getUseFpga();
if (use_gpu || use_fpga) {
auto wellContribs = WellContributions::create(accelerator_mode_, useWellConn);
bridge_->initWellContributions(*wellContribs, x.N() * x[0].N());
// the WellContributions can only be applied separately with CUDA or OpenCL, not with an FPGA or amgcl
#if HAVE_CUDA || HAVE_OPENCL
if (!useWellConn) {
getContribs(*wellContribs);
}
#endif
if (numJacobiBlocks_ > 1) {
this->copyMatToBlockJac(matrix, *blockJacobiForGPUILU0_);
// Const_cast needed since the CUDA stuff overwrites values for better matrix condition..
bridge_->solve_system(&matrix, blockJacobiForGPUILU0_.get(),
numJacobiBlocks_, rhs, *wellContribs, result);
}
else
bridge_->solve_system(&matrix, &matrix,
numJacobiBlocks_, rhs, *wellContribs, result);
if (result.converged) {
// get result vector x from non-Dune backend, iff solve was successful
bridge_->get_result(x);
return true;
} else {
// warn about CPU fallback
// BdaBridge might have disabled its BdaSolver for this simulation due to some error
// in that case the BdaBridge is disabled and flexibleSolver is always used
// or maybe the BdaSolver did not converge in time, then it will be used next linear solve
if (rank == 0) {
OpmLog::warning(bridge_->getAccleratorName() + " did not converge, now trying Dune to solve current linear system...");
}
}
}
return false;
}
template<class Matrix, class Vector>
template<class Grid>
void BdaSolverInfo<Matrix,Vector>::
blockJacobiAdjacency(const Grid& grid,
const std::vector<int>& cell_part,
size_t nonzeroes)
{
using size_type = typename Matrix::size_type;
using Iter = typename Matrix::CreateIterator;
size_type numCells = grid.size(0);
blockJacobiForGPUILU0_ = std::make_unique<Matrix>(numCells, numCells,
nonzeroes, Matrix::row_wise);
const auto& lid = grid.localIdSet();
const auto& gridView = grid.leafGridView();
auto elemIt = gridView.template begin<0>(); // should never overrun, since blockJacobiForGPUILU0_ is initialized with numCells rows
//Loop over cells
for (Iter row = blockJacobiForGPUILU0_->createbegin(); row != blockJacobiForGPUILU0_->createend(); ++elemIt, ++row)
{
const auto& elem = *elemIt;
size_type idx = lid.id(elem);
row.insert(idx);
// Add well non-zero connections
for (const auto wc : wellConnectionsGraph_[idx]) {
row.insert(wc);
}
int locPart = cell_part[idx];
//Add neighbor if it is on the same part
auto isend = gridView.iend(elem);
for (auto is = gridView.ibegin(elem); is!=isend; ++is)
{
//check if face has neighbor
if (is->neighbor())
{
size_type nid = lid.id(is->outside());
int nabPart = cell_part[nid];
if (locPart == nabPart) {
row.insert(nid);
}
}
}
}
}
template<class Matrix, class Vector>
void BdaSolverInfo<Matrix,Vector>::
copyMatToBlockJac(const Matrix& mat, Matrix& blockJac)
{
auto rbegin = blockJac.begin();
auto rend = blockJac.end();
auto outerRow = mat.begin();
for (auto row = rbegin; row != rend; ++row, ++outerRow) {
auto outerCol = (*outerRow).begin();
for (auto col = (*row).begin(); col != (*row).end(); ++col) {
// outerRow is guaranteed to have all column entries that row has!
while(outerCol.index() < col.index()) ++outerCol;
assert(outerCol.index() == col.index());
*col = *outerCol; // copy nonzero block
}
}
}
#endif
template<int Dim>
using BM = Dune::BCRSMatrix<MatrixBlock<double,Dim,Dim>>;
template<int Dim>
using BV = Dune::BlockVector<Dune::FieldVector<double,Dim>>;
#if HAVE_MPI
using CommunicationType = Dune::OwnerOverlapCopyCommunication<int,int>;
#else
using CommunicationType = Dune::CollectiveCommunication<int>;
#endif
#define INSTANCE_FLEX(Dim) \
template void makeOverlapRowsInvalid<BM<Dim>>(BM<Dim>&, const std::vector<int>&); \
template struct FlexibleSolverInfo<BM<Dim>,BV<Dim>,CommunicationType>;
#if HAVE_CUDA || HAVE_OPENCL || HAVE_FPGA || HAVE_AMGCL
#define INSTANCE(Dim) \
template struct BdaSolverInfo<BM<Dim>,BV<Dim>>; \
template void BdaSolverInfo<BM<Dim>,BV<Dim>>:: \
prepare<Dune::CpGrid>(const Dune::CpGrid&, \
const Dune::CartesianIndexMapper<Dune::CpGrid>&, \
const std::vector<Well>&, \
const std::vector<int>&, \
const size_t, const bool); \
INSTANCE_FLEX(Dim)
#else
#define INSTANCE(Dim) \
INSTANCE_FLEX(Dim)
#endif
INSTANCE(1)
INSTANCE(2)
INSTANCE(3)
INSTANCE(4)
INSTANCE(5)
INSTANCE(6)
}
}

457
opm/simulators/linalg/ISTLSolverEbos.hpp
View File

@ -22,24 +22,37 @@
#ifndef OPM_ISTLSOLVER_EBOS_HEADER_INCLUDED
#define OPM_ISTLSOLVER_EBOS_HEADER_INCLUDED
#include <dune/istl/owneroverlapcopy.hh>
#include <ebos/eclbasevanguard.hh>
#include <opm/common/ErrorMacros.hpp>
#include <opm/models/discretization/common/fvbaseproperties.hh>
#include <opm/models/common/multiphasebaseproperties.hh>
#include <opm/models/utils/parametersystem.hh>
#include <opm/models/utils/propertysystem.hh>
#include <opm/simulators/flow/BlackoilModelParametersEbos.hpp>
#include <opm/simulators/linalg/ExtractParallelGridInformationToISTL.hpp>
#include <opm/simulators/linalg/FlexibleSolver.hpp>
#include <opm/simulators/linalg/FlowLinearSolverParameters.hpp>
#include <opm/simulators/linalg/matrixblock.hh>
#include <opm/simulators/linalg/ParallelIstlInformation.hpp>
#include <opm/simulators/linalg/ParallelOverlappingILU0.hpp>
#include <opm/simulators/linalg/istlsparsematrixadapter.hh>
#include <opm/simulators/linalg/PreconditionerWithUpdate.hpp>
#include <opm/simulators/linalg/WellOperators.hpp>
#include <opm/simulators/linalg/WriteSystemMatrixHelper.hpp>
#include <opm/simulators/linalg/findOverlapRowsAndColumns.hpp>
#include <opm/simulators/linalg/getQuasiImpesWeights.hpp>
#include <opm/simulators/linalg/setupPropertyTree.hpp>
#if HAVE_CUDA || HAVE_OPENCL || HAVE_FPGA || HAVE_AMGCL
#include <opm/simulators/linalg/bda/BdaBridge.hpp>
#include <opm/simulators/linalg/bda/WellContributions.hpp>
#endif
#include <any>
#include <cstddef>
#include <functional>
#include <memory>
#include <set>
#include <sstream>
#include <string>
#include <tuple>
#include <vector>
namespace Opm::Properties {
@ -60,10 +73,10 @@ struct SparseMatrixAdapter<TypeTag, TTag::FlowIstlSolver>
private:
using Scalar = GetPropType<TypeTag, Properties::Scalar>;
enum { numEq = getPropValue<TypeTag, Properties::NumEq>() };
typedef MatrixBlock<Scalar, numEq, numEq> Block;
using Block = MatrixBlock<Scalar, numEq, numEq>;
public:
typedef typename Linear::IstlSparseMatrixAdapter<Block> type;
using type = typename Linear::IstlSparseMatrixAdapter<Block>;
};
} // namespace Opm::Properties
@ -71,6 +84,110 @@ public:
namespace Opm
{
#if HAVE_CUDA || HAVE_OPENCL || HAVE_FPGA || HAVE_AMGCL
template<class Matrix, class Vector, int block_size> class BdaBridge;
class WellContributions;
#endif
namespace detail
{
template<class Matrix, class Vector, class Comm>
struct FlexibleSolverInfo
{
using AbstractSolverType = Dune::InverseOperator<Vector, Vector>;
using AbstractOperatorType = Dune::AssembledLinearOperator<Matrix, Vector, Vector>;
using AbstractPreconditionerType = Dune::PreconditionerWithUpdate<Vector, Vector>;
void create(const Matrix& matrix,
bool parallel,
const PropertyTree& prm,
size_t pressureIndex,
std::function<Vector()> trueFunc,
Comm& comm);
std::unique_ptr<AbstractSolverType> solver_;
std::unique_ptr<AbstractOperatorType> op_;
std::unique_ptr<LinearOperatorExtra<Vector,Vector>> wellOperator_;
AbstractPreconditionerType* pre_ = nullptr;
size_t interiorCellNum_ = 0;
};
#if HAVE_CUDA || HAVE_OPENCL || HAVE_FPGA || HAVE_AMGCL
template<class Matrix, class Vector>
struct BdaSolverInfo
{
using WellContribFunc = std::function<void(WellContributions&)>;
using Bridge = BdaBridge<Matrix,Vector,Matrix::block_type::rows>;
BdaSolverInfo(const std::string& accelerator_mode,
const std::string& fpga_bitstream,
const int linear_solver_verbosity,
const int maxit,
const double tolerance,
const int platformID,
const int deviceID,
const std::string& opencl_ilu_reorder,
const std::string& linsolver);
~BdaSolverInfo();
template<class Grid>
void prepare(const Grid& grid,
const Dune::CartesianIndexMapper<Grid>& cartMapper,
const std::vector<Well>& wellsForConn,
const std::vector<int>& cellPartition,
const size_t nonzeroes,
const bool useWellConn);
bool apply(Vector& rhs,
const bool useWellConn,
WellContribFunc getContribs,
const int rank,
Matrix& matrix,
Vector& x,
Dune::InverseOperatorResult& result);
int numJacobiBlocks_ = 0;
private:
/// Create sparsity pattern for block-Jacobi matrix based on partitioning of grid.
/// Do not initialize the values, that is done in copyMatToBlockJac()
template<class Grid>
void blockJacobiAdjacency(const Grid& grid,
const std::vector<int>& cell_part,
size_t nonzeroes);
void copyMatToBlockJac(const Matrix& mat, Matrix& blockJac);
std::unique_ptr<Bridge> bridge_;
std::string accelerator_mode_;
std::unique_ptr<Matrix> blockJacobiForGPUILU0_;
std::vector<std::set<int>> wellConnectionsGraph_;
};
#endif
#ifdef HAVE_MPI
/// Copy values in parallel.
void copyParValues(std::any& parallelInformation, size_t size,
Dune::OwnerOverlapCopyCommunication<int,int>& comm);
#endif
/// Zero out off-diagonal blocks on rows corresponding to overlap cells
/// Diagonal blocks on ovelap rows are set to diag(1.0).
template<class Matrix>
void makeOverlapRowsInvalid(Matrix& matrix,
const std::vector<int>& overlapRows);
/// Create sparsity pattern for block-Jacobi matrix based on partitioning of grid.
/// Do not initialize the values, that is done in copyMatToBlockJac()
template<class Matrix, class Grid>
std::unique_ptr<Matrix> blockJacobiAdjacency(const Grid& grid,
const std::vector<int>& cell_part,
size_t nonzeroes,
const std::vector<std::set<int>>& wellConnectionsGraph);
}
/// This class solves the fully implicit black-oil system by
/// solving the reduced system (after eliminating well variables)
/// as a block-structured matrix (one block for all cell variables) for a fixed
@ -96,15 +213,10 @@ namespace Opm
using ElementMapper = GetPropType<TypeTag, Properties::ElementMapper>;
constexpr static std::size_t pressureIndex = GetPropType<TypeTag, Properties::Indices>::pressureSwitchIdx;
#if HAVE_CUDA || HAVE_OPENCL || HAVE_FPGA || HAVE_AMGCL
static const unsigned int block_size = Matrix::block_type::rows;
std::unique_ptr<BdaBridge<Matrix, Vector, block_size>> bdaBridge;
#endif
#if HAVE_MPI
using CommunicationType = Dune::OwnerOverlapCopyCommunication<int,int>;
#else
using CommunicationType = Dune::CollectiveCommunication< int >;
using CommunicationType = Dune::CollectiveCommunication<int>;
#endif
public:
@ -151,7 +263,15 @@ namespace Opm
const int linear_solver_verbosity = parameters_.linear_solver_verbosity_;
std::string fpga_bitstream = EWOMS_GET_PARAM(TypeTag, std::string, FpgaBitstream);
std::string linsolver = EWOMS_GET_PARAM(TypeTag, std::string, Linsolver);
bdaBridge.reset(new BdaBridge<Matrix, Vector, block_size>(accelerator_mode, fpga_bitstream, linear_solver_verbosity, maxit, tolerance, platformID, deviceID, opencl_ilu_reorder, linsolver));
bdaBridge = std::make_unique<detail::BdaSolverInfo<Matrix,Vector>>(accelerator_mode,
fpga_bitstream,
linear_solver_verbosity,
maxit,
tolerance,
platformID,
deviceID,
opencl_ilu_reorder,
linsolver);
}
#else
if (EWOMS_GET_PARAM(TypeTag, std::string, AcceleratorMode) != "none") {
@ -181,7 +301,7 @@ namespace Opm
OPM_THROW_NOLOG(std::runtime_error, msg);
}
interiorCellNum_ = detail::numMatrixRowsToUseInSolver(simulator_.vanguard().grid(), true);
flexibleSolver_.interiorCellNum_ = detail::numMatrixRowsToUseInSolver(simulator_.vanguard().grid(), true);
// Print parameters to PRT/DBG logs.
if (on_io_rank) {
@ -193,19 +313,17 @@ namespace Opm
}
// nothing to clean here
void eraseMatrix() {
void eraseMatrix()
{
}
void prepare(const SparseMatrixAdapter& M, Vector& b)
{
static bool firstcall = true;
#if HAVE_MPI
if (firstcall && parallelInformation_.type() == typeid(ParallelISTLInformation)) {
// Parallel case.
const ParallelISTLInformation* parinfo = std::any_cast<ParallelISTLInformation>(&parallelInformation_);
assert(parinfo);
if (firstcall) {
const size_t size = M.istlMatrix().N();
parinfo->copyValuesTo(comm_->indexSet(), comm_->remoteIndices(), size, 1);
detail::copyParValues(parallelInformation_, size, *comm_);
}
#endif
@ -216,22 +334,16 @@ namespace Opm
// Outch! We need to be able to scale the linear system! Hence const_cast
matrix_ = const_cast<Matrix*>(&M.istlMatrix());
useWellConn_ = EWOMS_GET_PARAM(TypeTag, bool, MatrixAddWellContributions);
// setup sparsity pattern for jacobi matrix for preconditioner (only used for openclSolver)
#if HAVE_OPENCL
this->numJacobiBlocks_ = EWOMS_GET_PARAM(TypeTag, int, NumJacobiBlocks);
#else
this->numJacobiBlocks_ = 0;
bdaBridge->numJacobiBlocks_ = EWOMS_GET_PARAM(TypeTag, int, NumJacobiBlocks);
bdaBridge->prepare(simulator_.vanguard().grid(),
simulator_.vanguard().cartesianIndexMapper(),
simulator_.vanguard().schedule().getWellsatEnd(),
simulator_.vanguard().cellPartition(),
getMatrix().nonzeroes(), useWellConn_);
#endif
useWellConn_ = EWOMS_GET_PARAM(TypeTag, bool, MatrixAddWellContributions);
if (numJacobiBlocks_ > 1) {
const auto wellsForConn = simulator_.vanguard().schedule().getWellsatEnd();
const auto& cartMapper = simulator_.vanguard().cartesianIndexMapper();
detail::setWellConnections(simulator_.vanguard().grid(), cartMapper, wellsForConn, useWellConn_,
wellConnectionsGraph_, numJacobiBlocks_);
std::cout << "Create block-Jacobi pattern" << std::endl;
blockJacobiAdjacency();
}
} else {
// Pointers should not change
if ( &(M.istlMatrix()) != matrix_ ) {
@ -242,26 +354,30 @@ namespace Opm
rhs_ = &b;
if (isParallel() && prm_.get<std::string>("preconditioner.type") != "ParOverILU0") {
makeOverlapRowsInvalid(getMatrix());
detail::makeOverlapRowsInvalid(getMatrix(), overlapRows_);
}
prepareFlexibleSolver();
firstcall = false;
}
void setResidual(Vector& /* b */) {
void setResidual(Vector& /* b */)
{
// rhs_ = &b; // Must be handled in prepare() instead.
}
void getResidual(Vector& b) const {
void getResidual(Vector& b) const
{
b = *rhs_;
}
void setMatrix(const SparseMatrixAdapter& /* M */) {
void setMatrix(const SparseMatrixAdapter& /* M */)
{
// matrix_ = &M.istlMatrix(); // Must be handled in prepare() instead.
}
bool solve(Vector& x) {
bool solve(Vector& x)
{
calls_ += 1;
// Write linear system if asked for.
const int verbosity = prm_.get<int>("verbosity", 0);
@ -275,54 +391,23 @@ namespace Opm
// Solve system.
Dune::InverseOperatorResult result;
bool accelerator_was_used = false;
// Use GPU if: available, chosen by user, and successful.
// Use FPGA if: support compiled, chosen by user, and successful.
#if HAVE_CUDA || HAVE_OPENCL || HAVE_FPGA || HAVE_AMGCL
bool use_gpu = bdaBridge->getUseGpu();
bool use_fpga = bdaBridge->getUseFpga();
if (use_gpu || use_fpga) {
const std::string accelerator_mode = EWOMS_GET_PARAM(TypeTag, std::string, AcceleratorMode);
auto wellContribs = WellContributions::create(accelerator_mode, useWellConn_);
bdaBridge->initWellContributions(*wellContribs, x.N() * x[0].N());
// the WellContributions can only be applied separately with CUDA or OpenCL, not with an FPGA or amgcl
#if HAVE_CUDA || HAVE_OPENCL
if (!useWellConn_) {
simulator_.problem().wellModel().getWellContributions(*wellContribs);
}
std::function<void(WellContributions&)> getContribs =
[this](WellContributions& w)
{
this->simulator_.problem().wellModel().getWellContributions(w);
};
if (!bdaBridge->apply(*rhs_, useWellConn_, getContribs,
simulator_.gridView().comm().rank(),
const_cast<Matrix&>(this->getMatrix()),
x, result))
#endif
if (numJacobiBlocks_ > 1) {
copyMatToBlockJac(getMatrix(), *blockJacobiForGPUILU0_);
// Const_cast needed since the CUDA stuff overwrites values for better matrix condition..
bdaBridge->solve_system(const_cast<Matrix*>(&getMatrix()), &*blockJacobiForGPUILU0_,
numJacobiBlocks_, *rhs_, *wellContribs, result);
}
else
bdaBridge->solve_system(const_cast<Matrix*>(&getMatrix()), const_cast<Matrix*>(&getMatrix()),
numJacobiBlocks_, *rhs_, *wellContribs, result);
if (result.converged) {
// get result vector x from non-Dune backend, iff solve was successful
bdaBridge->get_result(x);
accelerator_was_used = true;
} else {
// warn about CPU fallback
// BdaBridge might have disabled its BdaSolver for this simulation due to some error
// in that case the BdaBridge is disabled and flexibleSolver is always used
// or maybe the BdaSolver did not converge in time, then it will be used next linear solve
if (simulator_.gridView().comm().rank() == 0) {
OpmLog::warning(bdaBridge->getAccleratorName() + " did not converge, now trying Dune to solve current linear system...");
}
}
}
#endif
// Otherwise, use flexible istl solver.
if (!accelerator_was_used) {
assert(flexibleSolver_);
flexibleSolver_->apply(x, *rhs_, result);
{
assert(flexibleSolver_.solver_);
flexibleSolver_.solver_->apply(x, *rhs_, result);
}
// Check convergence, iterations etc.
@ -346,10 +431,7 @@ namespace Opm
protected:
#if HAVE_MPI
typedef Dune::OwnerOverlapCopyCommunication<int, int> Comm;
// 3x3 matrix block inversion was unstable from at least 2.3 until and
// including 2.5.0
typedef ParallelOverlappingILU0<Matrix,Vector,Vector,Comm> ParPreconditioner;
using Comm = Dune::OwnerOverlapCopyCommunication<int, int>;
#endif
void checkConvergence( const Dune::InverseOperatorResult& result ) const
@ -377,54 +459,28 @@ namespace Opm
void prepareFlexibleSolver()
{
std::function<Vector()> weightsCalculator = getWeightsCalculator();
if (shouldCreateSolver()) {
if (isParallel()) {
#if HAVE_MPI
if (useWellConn_) {
using ParOperatorType = Dune::OverlappingSchwarzOperator<Matrix, Vector, Vector, Comm>;
auto op = std::make_unique<ParOperatorType>(getMatrix(), *comm_);
using FlexibleSolverType = Dune::FlexibleSolver<ParOperatorType>;
auto sol = std::make_unique<FlexibleSolverType>(*op, *comm_, prm_, weightsCalculator, pressureIndex);
preconditionerForFlexibleSolver_ = &(sol->preconditioner());
linearOperatorForFlexibleSolver_ = std::move(op);
flexibleSolver_ = std::move(sol);
} else {
using ParOperatorType = WellModelGhostLastMatrixAdapter<Matrix, Vector, Vector, true>;
wellOperator_ = std::make_unique<WellModelOperator>(simulator_.problem().wellModel());
auto op = std::make_unique<ParOperatorType>(getMatrix(), *wellOperator_, interiorCellNum_);
using FlexibleSolverType = Dune::FlexibleSolver<ParOperatorType>;
auto sol = std::make_unique<FlexibleSolverType>(*op, *comm_, prm_, weightsCalculator, pressureIndex);
preconditionerForFlexibleSolver_ = &(sol->preconditioner());
linearOperatorForFlexibleSolver_ = std::move(op);
flexibleSolver_ = std::move(sol);
}
#endif
} else {
if (useWellConn_) {
using SeqOperatorType = Dune::MatrixAdapter<Matrix, Vector, Vector>;
auto op = std::make_unique<SeqOperatorType>(getMatrix());
using FlexibleSolverType = Dune::FlexibleSolver<SeqOperatorType>;
auto sol = std::make_unique<FlexibleSolverType>(*op, prm_, weightsCalculator, pressureIndex);
preconditionerForFlexibleSolver_ = &(sol->preconditioner());
linearOperatorForFlexibleSolver_ = std::move(op);
flexibleSolver_ = std::move(sol);
} else {
using SeqOperatorType = WellModelMatrixAdapter<Matrix, Vector, Vector, false>;
wellOperator_ = std::make_unique<WellModelOperator>(simulator_.problem().wellModel());
auto op = std::make_unique<SeqOperatorType>(getMatrix(), *wellOperator_);
using FlexibleSolverType = Dune::FlexibleSolver<SeqOperatorType>;
auto sol = std::make_unique<FlexibleSolverType>(*op, prm_, weightsCalculator, pressureIndex);
preconditionerForFlexibleSolver_ = &(sol->preconditioner());
linearOperatorForFlexibleSolver_ = std::move(op);
flexibleSolver_ = std::move(sol);
}
std::function<Vector()> trueFunc =
[this]
{
return this->getTrueImpesWeights(pressureIndex);
};
if (!useWellConn_) {
auto wellOp = std::make_unique<WellModelOperator>(simulator_.problem().wellModel());
flexibleSolver_.wellOperator_ = std::move(wellOp);
}
flexibleSolver_.create(getMatrix(),
isParallel(),
prm_,
pressureIndex,
trueFunc,
*comm_);
}
else
{
preconditionerForFlexibleSolver_->update();
flexibleSolver_.pre_->update();
}
}
@ -435,7 +491,7 @@ namespace Opm
{
// Decide if we should recreate the solver or just do
// a minimal preconditioner update.
if (!flexibleSolver_) {
if (!flexibleSolver_.solver_) {
return true;
}
if (this->parameters_.cpr_reuse_setup_ == 0) {
@ -476,41 +532,6 @@ namespace Opm
}
/// Return an appropriate weight function if a cpr preconditioner is asked for.
std::function<Vector()> getWeightsCalculator() const
{
std::function<Vector()> weightsCalculator;
using namespace std::string_literals;
auto preconditionerType = prm_.get("preconditioner.type"s, "cpr"s);
if (preconditionerType == "cpr" || preconditionerType == "cprt"
|| preconditionerType == "cprw" || preconditionerType == "cprwt") {
const bool transpose = preconditionerType == "cprt" || preconditionerType == "cprwt";
const auto weightsType = prm_.get("preconditioner.weight_type"s, "quasiimpes"s);
if (weightsType == "quasiimpes") {
// weights will be created as default in the solver
// assignment p = pressureIndex prevent compiler warning about
// capturing variable with non-automatic storage duration
weightsCalculator = [this, transpose, p = pressureIndex]() {
return Amg::getQuasiImpesWeights<Matrix, Vector>(this->getMatrix(), p, transpose);
};
} else if (weightsType == "trueimpes") {
// assignment p = pressureIndex prevent compiler warning about
// capturing variable with non-automatic storage duration
weightsCalculator = [this, p = pressureIndex]() {
return this->getTrueImpesWeights(p);
};
} else {
OPM_THROW(std::invalid_argument,
"Weights type " << weightsType << "not implemented for cpr."
<< " Please use quasiimpes or trueimpes.");
}
}
return weightsCalculator;
}
// Weights to make approximate pressure equations.
// Calculated from the storage terms (only) of the
// conservation equations, ignoring all other terms.
@ -518,98 +539,13 @@ namespace Opm
{
Vector weights(rhs_->size());
ElementContext elemCtx(simulator_);
Amg::getTrueImpesWeights(pressureVarIndex, weights, simulator_.vanguard().gridView(),
Amg::getTrueImpesWeights(pressureVarIndex, weights,
simulator_.vanguard().gridView(),
elemCtx, simulator_.model(),
ThreadManager::threadId());
return weights;
}
/// Zero out off-diagonal blocks on rows corresponding to overlap cells
/// Diagonal blocks on ovelap rows are set to diag(1.0).
void makeOverlapRowsInvalid(Matrix& matrix) const
{
//value to set on diagonal
const int numEq = Matrix::block_type::rows;
typename Matrix::block_type diag_block(0.0);
for (int eq = 0; eq < numEq; ++eq)
diag_block[eq][eq] = 1.0;
//loop over precalculated overlap rows and columns
for (auto row = overlapRows_.begin(); row != overlapRows_.end(); row++ )
{
int lcell = *row;
// Zero out row.
matrix[lcell] = 0.0;
//diagonal block set to diag(1.0).
matrix[lcell][lcell] = diag_block;
}
}
/// Create sparsity pattern for block-Jacobi matrix based on partitioning of grid.
/// Do not initialize the values, that is done in copyMatToBlockJac()
void blockJacobiAdjacency()
{
const auto& grid = simulator_.vanguard().grid();
std::vector<int> cell_part = simulator_.vanguard().cellPartition();
typedef typename Matrix::size_type size_type;
typedef typename Matrix::CreateIterator Iter;
size_type numCells = grid.size( 0 );
blockJacobiForGPUILU0_.reset(new Matrix(numCells, numCells, getMatrix().nonzeroes(), Matrix::row_wise));
const auto& lid = grid.localIdSet();
const auto& gridView = grid.leafGridView();
auto elemIt = gridView.template begin<0>(); // should never overrun, since blockJacobiForGPUILU0_ is initialized with numCells rows
//Loop over cells
for (Iter row = blockJacobiForGPUILU0_->createbegin(); row != blockJacobiForGPUILU0_->createend(); ++elemIt, ++row)
{
const auto& elem = *elemIt;
size_type idx = lid.id(elem);
row.insert(idx);
// Add well non-zero connections
for (auto wc = wellConnectionsGraph_[idx].begin(); wc!=wellConnectionsGraph_[idx].end(); ++wc) {
row.insert(*wc);
}
int locPart = cell_part[idx];
//Add neighbor if it is on the same part
auto isend = gridView.iend(elem);
for (auto is = gridView.ibegin(elem); is!=isend; ++is)
{
//check if face has neighbor
if (is->neighbor())
{
size_type nid = lid.id(is->outside());
int nabPart = cell_part[nid];
if (locPart == nabPart) {
row.insert(nid);
}
}
}
}
}
void copyMatToBlockJac(const Matrix& mat, Matrix& blockJac)
{
auto rbegin = blockJac.begin();
auto rend = blockJac.end();
auto outerRow = mat.begin();
for (auto row = rbegin; row != rend; ++row, ++outerRow) {
auto outerCol = (*outerRow).begin();
for (auto col = (*row).begin(); col != (*row).end(); ++col) {
// outerRow is guaranteed to have all column entries that row has!
while(outerCol.index() < col.index()) ++outerCol;
assert(outerCol.index() == col.index());
*col = *outerCol; // copy nonzero block
}
}
}
Matrix& getMatrix()
{
return *matrix_;
@ -630,23 +566,18 @@ namespace Opm
Matrix* matrix_;
Vector *rhs_;
std::unique_ptr<Matrix> blockJacobiForGPUILU0_;
#if HAVE_CUDA || HAVE_OPENCL || HAVE_FPGA || HAVE_AMGCL
std::unique_ptr<detail::BdaSolverInfo<Matrix, Vector>> bdaBridge;
#endif
std::unique_ptr<AbstractSolverType> flexibleSolver_;
std::unique_ptr<AbstractOperatorType> linearOperatorForFlexibleSolver_;
AbstractPreconditionerType* preconditionerForFlexibleSolver_;
std::unique_ptr<WellModelAsLinearOperator<WellModel, Vector, Vector>> wellOperator_;
detail::FlexibleSolverInfo<Matrix,Vector,CommunicationType> flexibleSolver_;
std::vector<int> overlapRows_;
std::vector<int> interiorRows_;
std::vector<std::set<int>> wellConnectionsGraph_;
bool useWellConn_;
size_t interiorCellNum_;
FlowLinearSolverParameters parameters_;
PropertyTree prm_;
bool scale_variables_;
int numJacobiBlocks_;
std::shared_ptr< CommunicationType > comm_;
}; // end ISTLSolver

3
tests/test_glift1.cpp
View File

@ -40,10 +40,7 @@
#include <opm/simulators/wells/GasLiftSingleWell.hpp>
#include <opm/simulators/wells/GasLiftSingleWellGeneric.hpp>
#include <opm/simulators/wells/GasLiftGroupInfo.hpp>
//#include <opm/simulators/flow/SimulatorFullyImplicitBlackoilEbos.hpp>
//#include <flow/flow_ebos_blackoil.hpp>
#include <opm/simulators/wells/WellState.hpp>
/// #include <opm/simulators/flow/Main.hpp>
#if HAVE_DUNE_FEM
#include <dune/fem/misc/mpimanager.hh>