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

Changes to make CPR able to use general indexing.

Browse Source
This commit is contained in:
hnil 2019-03-15 11:27:38 +01:00 committed by Atgeirr Flø Rasmussen
parent 9af17dcc8c
commit 31f4c01938
4 changed files with 78 additions and 57 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

67
opm/autodiff/BlackoilAmg.hpp
View File

@ -42,6 +42,27 @@ class UnSymmetricCriterion;
}
}
namespace Opm
{
namespace Amg
{
template<int Row, int Column>
class Element
{
public:
enum { /* @brief We preserve the sign.*/
is_sign_preserving = true
};
template<class M>
typename M::field_type operator()(const M& m) const
{
return m[Row][Column];
}
};
} // namespace Amg
} // namespace Opm
namespace Dune
{
@ -99,12 +120,12 @@ Dune::OverlappingSchwarzOperator<M,X,Y,T> createOperator(const Dune::Overlapping
//! Sedimentary Basin Simulations, 2003.
//! \param op The operator that stems from the discretization.
//! \param comm The communication objecte describing the data distribution.
//! \param pressureIndex The index of the pressure in the matrix block
//! \param pressureEqnIndex The index of the pressure in the matrix block
//! \retun A pair of the scaled matrix and the associated operator-
template<class Operator, class Communication, class Vector>
std::tuple<std::unique_ptr<typename Operator::matrix_type>, Operator>
scaleMatrixDRS(const Operator& op, const Communication& comm,
std::size_t pressureIndex, const Vector& weights, const Opm::CPRParameter& param)
std::size_t pressureEqnIndex, const Vector& weights, const Opm::CPRParameter& param)
{
using Matrix = typename Operator::matrix_type;
using Block = typename Matrix::block_type;
@ -117,7 +138,7 @@ scaleMatrixDRS(const Operator& op, const Communication& comm,
const auto endj = (*i).end();
for (auto j = (*i).begin(); j != endj; ++j) {
Block& block = *j;
BlockVector& bvec = block[pressureIndex];
BlockVector& bvec = block[pressureEqnIndex];
// should introduce limits which also change the weights
block.mtv(bw, bvec);
}
@ -131,16 +152,16 @@ scaleMatrixDRS(const Operator& op, const Communication& comm,
//! See section 3.2.3 of Scheichl, Masson: Decoupling and Block Preconditioning for
//! Sedimentary Basin Simulations, 2003.
//! \param vector The vector to scale
//! \param pressureIndex The index of the pressure in the matrix block
//! \param pressureEqnIndex The index of the pressure in the matrix block
template<class Vector>
void scaleVectorDRS(Vector& vector, std::size_t pressureIndex, const Opm::CPRParameter& param, const Vector& weights)
void scaleVectorDRS(Vector& vector, std::size_t pressureEqnIndex, const Opm::CPRParameter& param, const Vector& weights)
{
using Block = typename Vector::block_type;
if (param.cpr_use_drs_) {
for (std::size_t j = 0; j < vector.size(); ++j) {
Block& block = vector[j];
const Block& bw = weights[j];
block[pressureIndex] = bw.dot(block);
block[pressureEqnIndex] = bw.dot(block);
}
}
}
@ -289,23 +310,23 @@ struct ScalarType<Dune::Amg::CoarsenCriterion<Dune::Amg::UnSymmetricCriterion<Du
using value = Dune::Amg::CoarsenCriterion<Dune::Amg::UnSymmetricCriterion<Dune::BCRSMatrix<typename ScalarType<B>::value>, Dune::Amg::FirstDiagonal> >;
};
template<class C, std::size_t COMPONENT_INDEX>
template<class C, std::size_t COMPONENT_INDEX, std::size_t VARIABLE_INDEX>
struct OneComponentCriterionType
{};
template<class B, class N, std::size_t COMPONENT_INDEX>
struct OneComponentCriterionType<Dune::Amg::CoarsenCriterion<Dune::Amg::SymmetricCriterion<Dune::BCRSMatrix<B>,N> >,COMPONENT_INDEX>
template<class B, class N, std::size_t COMPONENT_INDEX, std::size_t VARIABLE_INDEX>
struct OneComponentCriterionType<Dune::Amg::CoarsenCriterion<Dune::Amg::SymmetricCriterion<Dune::BCRSMatrix<B>,N> >, COMPONENT_INDEX, VARIABLE_INDEX>
{
using value = Dune::Amg::CoarsenCriterion<Dune::Amg::SymmetricCriterion<Dune::BCRSMatrix<B>, Dune::Amg::Diagonal<COMPONENT_INDEX> > >;
using value = Dune::Amg::CoarsenCriterion<Dune::Amg::SymmetricCriterion<Dune::BCRSMatrix<B>, Opm::Amg::Element<COMPONENT_INDEX, VARIABLE_INDEX> > >;
};
template<class B, class N, std::size_t COMPONENT_INDEX>
struct OneComponentCriterionType<Dune::Amg::CoarsenCriterion<Dune::Amg::UnSymmetricCriterion<Dune::BCRSMatrix<B>,N> >,COMPONENT_INDEX>
template<class B, class N, std::size_t COMPONENT_INDEX, std::size_t VARIABLE_INDEX>
struct OneComponentCriterionType<Dune::Amg::CoarsenCriterion<Dune::Amg::UnSymmetricCriterion<Dune::BCRSMatrix<B>,N> >, COMPONENT_INDEX, VARIABLE_INDEX>
{
using value = Dune::Amg::CoarsenCriterion<Dune::Amg::UnSymmetricCriterion<Dune::BCRSMatrix<B>, Dune::Amg::Diagonal<COMPONENT_INDEX> > >;
using value = Dune::Amg::CoarsenCriterion<Dune::Amg::UnSymmetricCriterion<Dune::BCRSMatrix<B>, Opm::Amg::Element<COMPONENT_INDEX, VARIABLE_INDEX> > >;
};
template<class Operator, class Criterion, class Communication, std::size_t COMPONENT_INDEX>
template<class Operator, class Criterion, class Communication, std::size_t COMPONENT_INDEX, std::size_t VARIABLE_INDEX>
class OneComponentAggregationLevelTransferPolicy;
@ -674,7 +695,7 @@ void buildCoarseSparseMatrix(M& coarseMatrix, G& fineGraph, const V& visitedMap,
* @tparam Criterion The criterion that describes the aggregation procedure.
* @tparam Communication The class that describes the communication pattern.
*/
template<class Operator, class Criterion, class Communication, std::size_t COMPONENT_INDEX>
template<class Operator, class Criterion, class Communication, std::size_t COMPONENT_INDEX, std::size_t VARIABLE_INDEX>
class OneComponentAggregationLevelTransferPolicy
: public Dune::Amg::LevelTransferPolicy<Operator, typename Detail::ScalarType<Operator>::value>
{
@ -785,7 +806,7 @@ public:
for ( auto col = row.begin(), cend = row.end(); col != cend; ++col, ++coarseCol )
{
assert( col.index() == coarseCol.index() );
*coarseCol = (*col)[COMPONENT_INDEX][COMPONENT_INDEX];
*coarseCol = (*col)[COMPONENT_INDEX][VARIABLE_INDEX];
}
++coarseRow;
}
@ -815,7 +836,7 @@ public:
const auto& j = (*aggregatesMap_)[entry.index()];
if ( j != AggregatesMap::ISOLATED )
{
(*coarseLevelMatrix_)[i][j] += (*entry)[COMPONENT_INDEX][COMPONENT_INDEX];
(*coarseLevelMatrix_)[i][j] += (*entry)[COMPONENT_INDEX][VARIABLE_INDEX];
}
}
}
@ -911,9 +932,10 @@ private:
* \tparam C The type of coarsening criterion to use.
* \tparam P The type of the class describing the parallelization.
* \tparam COMPONENT_INDEX The index of the component to use for coarsening (usually the pressure).
* \tparam VARIABLE_INDEX The index of the variable to use for coarsening (usually the pressure).
*/
template<typename O, typename S, typename C,
typename P, std::size_t COMPONENT_INDEX>
typename P, std::size_t COMPONENT_INDEX, std::size_t VARIABLE_INDEX>
class BlackoilAmg
: public Dune::Preconditioner<typename O::domain_type, typename O::range_type>
{
@ -934,13 +956,14 @@ protected:
using CoarseOperator = typename Detail::ScalarType<Operator>::value;
using CoarseSmoother = typename Detail::ScalarType<Smoother>::value;
using FineCriterion =
typename Detail::OneComponentCriterionType<Criterion,COMPONENT_INDEX>::value;
typename Detail::OneComponentCriterionType<Criterion, COMPONENT_INDEX, VARIABLE_INDEX>::value;
using CoarseCriterion = typename Detail::ScalarType<Criterion>::value;
using LevelTransferPolicy =
OneComponentAggregationLevelTransferPolicy<Operator,
FineCriterion,
Communication,
COMPONENT_INDEX>;
COMPONENT_INDEX,
VARIABLE_INDEX>;
using CoarseSolverPolicy =
Detail::OneStepAMGCoarseSolverPolicy<CoarseOperator,
CoarseSmoother,
@ -1028,7 +1051,7 @@ namespace ISTLUtility
/// \tparam C The type of the coarsening criterion to use.
/// \tparam index The pressure index.
////
template<class M, class X, class Y, class P, class C, std::size_t index>
template<class M, class X, class Y, class P, class C, std::size_t pressureEqnIndex, std::size_t pressureVarIndex>
struct BlackoilAmgSelector
{
using Criterion = C;
@ -1036,7 +1059,7 @@ struct BlackoilAmgSelector
using ParallelInformation = typename Selector::ParallelInformation;
using Operator = typename Selector::Operator;
using Smoother = typename Selector::EllipticPreconditioner;
using AMG = BlackoilAmg<Operator,Smoother,Criterion,ParallelInformation,index>;
using AMG = BlackoilAmg<Operator, Smoother, Criterion, ParallelInformation, pressureEqnIndex, pressureVarIndex>;
};
} // end namespace ISTLUtility
} // end namespace Opm

18
opm/autodiff/CPRPreconditioner.hpp
View File

@ -55,9 +55,15 @@ namespace Opm
{
template<typename O, typename S, typename C,
typename P, std::size_t COMPONENT_INDEX>
typename P, std::size_t COMPONENT_INDEX, std::size_t VARIABLE_INDEX>
class BlackoilAmg;
namespace Amg
{
template<int Row, int Column>
class Element;
}
namespace ISTLUtility
{
@ -183,14 +189,14 @@ createEllipticPreconditionerPointer(const M& Ae, double relax,
return EllipticPreconditionerPointer(new ParallelPreconditioner(Ae, comm, relax, milu));
}
template < class C, class Op, class P, class S, std::size_t index, class Vector>
template < class C, class Op, class P, class S, std::size_t PressureEqnIndex, std::size_t PressureVarIndex, class Vector>
inline void
createAMGPreconditionerPointer(Op& opA, const double relax, const P& comm,
std::unique_ptr< BlackoilAmg<Op,S,C,P,index> >& amgPtr,
std::unique_ptr< BlackoilAmg<Op,S,C,P,PressureEqnIndex,PressureVarIndex> >& amgPtr,
const CPRParameter& params,
const Vector& weights)
{
using AMG = BlackoilAmg<Op,S,C,P,index>;
using AMG = BlackoilAmg<Op,S,C,P,PressureEqnIndex,PressureVarIndex>;
const int verbosity = ( params.cpr_solver_verbose_ && comm.communicator().rank() == 0 ) ? 1 : 0;
// TODO: revise choice of parameters
@ -242,7 +248,7 @@ createAMGPreconditionerPointer(Op& opA, const double relax, const MILU_VARIANT m
/// \param relax The relaxation parameter for ILU0.
/// \param comm The object describing the parallelization information and communication.
// \param amgPtr The unique_ptr to be filled (return)
template < int pressureIndex=0, class Op, class P, class AMG >
template < int PressureEqnIndex, int PressureVarIndex, class Op, class P, class AMG >
inline void
createAMGPreconditionerPointer( Op& opA, const double relax, const MILU_VARIANT milu, const P& comm, std::unique_ptr< AMG >& amgPtr )
{
@ -250,7 +256,7 @@ createAMGPreconditionerPointer( Op& opA, const double relax, const MILU_VARIANT
typedef typename Op::matrix_type M;
// The coupling metric used in the AMG
typedef Dune::Amg::Diagonal<pressureIndex> CouplingMetric;
typedef Opm::Amg::Element<PressureEqnIndex, PressureVarIndex> CouplingMetric;
// The coupling criterion used in the AMG
typedef Dune::Amg::SymmetricCriterion<M, CouplingMetric> CritBase;

48
opm/autodiff/ISTLSolverEbos.hpp
View File

@ -1,5 +1,6 @@
/*
Copyright 2016 IRIS AS
Copyright 2019 Equinor ASA
This file is part of the Open Porous Media project (OPM).
@ -31,6 +32,7 @@
#include <opm/common/Exceptions.hpp>
#include <opm/core/linalg/ParallelIstlInformation.hpp>
#include <opm/common/utility/platform_dependent/disable_warnings.h>
#include <opm/material/fluidsystems/BlackOilDefaultIndexTraits.hpp>
#include <ewoms/common/parametersystem.hh>
#include <ewoms/common/propertysystem.hh>
@ -162,11 +164,6 @@ protected:
/// solving the reduced system (after eliminating well variables)
/// as a block-structured matrix (one block for all cell variables) for a fixed
/// number of cell variables np .
/// \tparam MatrixBlockType The type of the matrix block used.
/// \tparam VectorBlockType The type of the vector block used.
/// \tparam pressureIndex The index of the pressure component in the vector
/// vector block. It is used to guide the AMG coarsening.
/// Default is zero.
template <class TypeTag>
class ISTLSolverEbos
{
@ -185,7 +182,9 @@ protected:
typedef typename GET_PROP_TYPE(TypeTag, ThreadManager) ThreadManager;
typedef typename GridView::template Codim<0>::Entity Element;
typedef typename GET_PROP_TYPE(TypeTag, ElementContext) ElementContext;
enum { pressureIndex = Indices::pressureSwitchIdx };
// Due to miscibility oil <-> gas the water eqn is the one we can replace with a pressure equation.
enum { pressureEqnIndex = BlackOilDefaultIndexTraits::waterCompIdx };
enum { pressureVarIndex = Indices::pressureSwitchIdx };
static const int numEq = Indices::numEq;
public:
@ -236,7 +235,7 @@ protected:
weights_ = getSimpleWeights(bvec);
} else if (parameters_.system_strategy_ == "original") {
BlockVector bvec(0.0);
bvec[pressureIndex] = 1;
bvec[pressureEqnIndex] = 1;
weights_ = getSimpleWeights(bvec);
} else {
form_cpr = false;
@ -370,11 +369,11 @@ protected:
if ( parameters_.use_cpr_ )
{
using Matrix = typename MatrixOperator::matrix_type;
using CouplingMetric = Dune::Amg::Diagonal<pressureIndex>;
using CouplingMetric = Opm::Amg::Element<pressureEqnIndex, pressureVarIndex>;
using CritBase = Dune::Amg::SymmetricCriterion<Matrix, CouplingMetric>;
using Criterion = Dune::Amg::CoarsenCriterion<CritBase>;
using AMG = typename ISTLUtility
::BlackoilAmgSelector< Matrix, Vector, Vector,POrComm, Criterion, pressureIndex >::AMG;
::BlackoilAmgSelector< Matrix, Vector, Vector,POrComm, Criterion, pressureEqnIndex, pressureVarIndex >::AMG;
std::unique_ptr< AMG > amg;
// Construct preconditioner.
@ -458,7 +457,7 @@ protected:
void
constructAMGPrecond(LinearOperator& /* linearOperator */, const POrComm& comm, std::unique_ptr< AMG >& amg, std::unique_ptr< MatrixOperator >& opA, const double relax, const MILU_VARIANT milu) const
{
ISTLUtility::template createAMGPreconditionerPointer<pressureIndex>( *opA, relax, milu, comm, amg );
ISTLUtility::template createAMGPreconditionerPointer<pressureEqnIndex, pressureVarIndex>( *opA, relax, milu, comm, amg );
}
@ -636,7 +635,7 @@ protected:
{
Vector weights(rhs_->size());
BlockVector rhs(0.0);
rhs[pressureIndex] = 1.0;
rhs[pressureVarIndex] = 1.0;
int index = 0;
ElementContext elemCtx(simulator_);
const auto& vanguard = simulator_.vanguard();
@ -657,7 +656,7 @@ protected:
for (int ii = 0; ii < numEq; ++ii) {
for (int jj = 0; jj < numEq; ++jj) {
block[ii][jj] = storage[ii].derivative(jj)/storage_scale;
if (jj == 0) {
if (jj == pressureVarIndex) {
block[ii][jj] *= pressure_scale;
}
}
@ -720,7 +719,7 @@ protected:
Matrix& A = *matrix_;
Vector weights(rhs_->size());
BlockVector rhs(0.0);
rhs[pressureIndex] = 1;
rhs[pressureVarIndex] = 1;
const auto endi = A.end();
for (auto i = A.begin(); i!=endi; ++i) {
const auto endj = (*i).end();
@ -761,27 +760,20 @@ protected:
const auto endj = (*i).end();
for (auto j = (*i).begin(); j != endj; ++j) {
// assume it is something on all rows
// the blew logic depend on pressureIndex=0
Block& block = (*j);
for ( std::size_t ii = 0; ii < block.rows; ii++ ) {
if ( ii == 0 ) {
for (std::size_t jj = 0; jj < block.cols; jj++) {
block[0][jj] *= bweights[ii];
}
} else {
for (std::size_t jj = 0; jj < block.cols; jj++) {
block[0][jj] += bweights[ii]*block[ii][jj];
}
BlockVector neweq(0.0);
for (std::size_t ii = 0; ii < block.rows; ii++) {
for (std::size_t jj = 0; jj < block.cols; jj++) {
neweq[jj] += bweights[ii]*block[ii][jj];
}
}
block[pressureEqnIndex] = neweq;
}
BlockVector newrhs(0.0);
for (std::size_t ii = 0; ii < brhs.size(); ii++) {
if ( ii == 0 ){
brhs[0] *= bweights[ii];
} else {
brhs[0] += bweights[ii]*brhs[ii];
}
newrhs += bweights[ii]*brhs[ii];
}
brhs = newrhs;
}
}

2
tests/test_blackoil_amg.cpp
View File

@ -303,7 +303,7 @@ void runBlackoilAmgLaplace()
smootherArgs.iterations = 1;
Opm::CPRParameter param;
Opm::BlackoilAmg<Operator,ParSmoother,Criterion,Communication,0> amg(param,
Opm::BlackoilAmg<Operator,ParSmoother,Criterion,Communication,0,0> amg(param,
{},
fop, criterion,
smootherArgs,