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

Removes stale CPR headers that have been superseeded.

Browse Source 
No need to drag them along and confuse people.
This commit is contained in:
Markus Blatt 2020-04-15 21:38:28 +02:00
parent 83d30547f6
commit 057a0ceeeb
3 changed files with 0 additions and 478 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

2
CMakeLists_files.cmake
View File

@ -146,7 +146,6 @@ list (APPEND PUBLIC_HEADER_FILES
opm/simulators/linalg/bda/cusparseSolverBackend.hpp
opm/simulators/linalg/bda/WellContributions.hpp
opm/simulators/linalg/BlackoilAmg.hpp
opm/simulators/linalg/BlackoilAmgCpr.hpp
opm/simulators/linalg/amgcpr.hh
opm/simulators/linalg/twolevelmethodcpr.hh
opm/simulators/linalg/CPRPreconditioner.hpp
@ -155,7 +154,6 @@ list (APPEND PUBLIC_HEADER_FILES
opm/simulators/linalg/FlowLinearSolverParameters.hpp
opm/simulators/linalg/GraphColoring.hpp
opm/simulators/linalg/ISTLSolverEbos.hpp
opm/simulators/linalg/ISTLSolverEbosCpr.hpp
opm/simulators/linalg/ISTLSolverEbosFlexible.hpp
opm/simulators/linalg/MatrixBlock.hpp
opm/simulators/linalg/OwningBlockPreconditioner.hpp

173
opm/simulators/linalg/BlackoilAmgCpr.hpp
View File

@ -1,173 +0,0 @@
/*
Copyright 2017 Dr. Blatt - HPC-Simulation-Software & Services
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_AMGCPR_HEADER_INCLUDED
#define OPM_AMGCPR_HEADER_INCLUDED
#include <opm/simulators/linalg/twolevelmethodcpr.hh>
#include <opm/simulators/linalg/matrixblock.hh>
#include <opm/simulators/linalg/ParallelOverlappingILU0.hpp>
#include <opm/simulators/linalg/FlowLinearSolverParameters.hpp>
#include <opm/simulators/linalg/CPRPreconditioner.hpp>
#include <opm/simulators/linalg/amgcpr.hh>
#include <dune/istl/paamg/twolevelmethod.hh>
#include <dune/istl/paamg/aggregates.hh>
#include <dune/istl/bvector.hh>
#include <dune/istl/bcrsmatrix.hh>
#include <dune/istl/preconditioners.hh>
#include <dune/istl/schwarz.hh>
#include <dune/istl/operators.hh>
#include <dune/istl/scalarproducts.hh>
#include <dune/common/fvector.hh>
#include <dune/common/fmatrix.hh>
namespace Opm
{
/**
* \brief An algebraic twolevel or multigrid approach for solving blackoil (supports CPR with and without AMG)
*
* This preconditioner first decouples the component used for coarsening using a simple scaling
* approach (e.g. Scheichl, Masson 2013,\see scaleMatrixDRS). Then it constructs the
* coarse level system. The coupling is defined by the weights corresponding to the element located at
* (COMPONENT_INDEX, VARIABLE_INDEX) in the block matrix. Then the coarse level system is constructed
* either by extracting these elements, or by applying aggregation to them directly. This coarse level
* can be solved either by AMG or by ILU. The preconditioner is configured using CPRParameter.
* \tparam O The type of the operator (encapsulating a BCRSMatrix).
* \tparam S The type of the smoother.
* \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 water).
* \tparam VARIABLE_INDEX The index of the variable to use for coarsening (usually pressure).
*/
template<typename O, typename S, typename SC, typename C,
typename P, std::size_t COMPONENT_INDEX, std::size_t VARIABLE_INDEX>
class BlackoilAmgCpr
: public Dune::Preconditioner<typename O::domain_type, typename O::range_type>
{
public:
/** \brief The type of the operator (encapsulating a BCRSMatrix). */
using Operator = O;
/** \brief The type of coarsening criterion to use. */
using Criterion = C;
/** \brief The type of the class describing the parallelization. */
using Communication = P;
/** \brief The type of the smoother. */
using Smoother = S;
/** \brief The type of the smoother arguments for construction. */
using SmootherArgs = typename Dune::Amg::SmootherTraits<Smoother>::Arguments;
protected:
using Matrix = typename Operator::matrix_type;
using CoarseOperator = typename Detail::ScalarType<Operator>::value;
using CoarseSmoother = typename Detail::ScalarType<SC>::value;
using FineCriterion =
typename Detail::OneComponentCriterionType<Criterion,COMPONENT_INDEX, VARIABLE_INDEX>::value;
using CoarseCriterion = typename Detail::ScalarType<Criterion>::value;
using LevelTransferPolicy =
OneComponentAggregationLevelTransferPolicy<Operator,
FineCriterion,
Communication,
COMPONENT_INDEX,
VARIABLE_INDEX>;
using CoarseSolverPolicy =
Detail::OneStepAMGCoarseSolverPolicy<CoarseOperator,
CoarseSmoother,
CoarseCriterion,
LevelTransferPolicy>;
using TwoLevelMethod =
Dune::Amg::TwoLevelMethodCpr<Operator,
CoarseSolverPolicy,
Smoother>;
public:
Dune::SolverCategory::Category category() const override
{
return std::is_same<Communication, Dune::Amg::SequentialInformation>::value ?
Dune::SolverCategory::sequential : Dune::SolverCategory::overlapping;
}
/**
* \brief Constructor.
* \param param The parameters used for configuring the solver.
* \param fineOperator The operator of the fine level.
* \param criterion The criterion describing the coarsening approach.
* \param smargs The arguments for constructing the smoother.
* \param comm The information about the parallelization.
*/
BlackoilAmgCpr(const CPRParameter& param,
const typename TwoLevelMethod::FineDomainType& weights,
const Operator& fineOperator, const Criterion& criterion,
const SmootherArgs& smargs, const Communication& comm)
: param_(param),
weights_(weights),
scaledMatrix_(Detail::scaleMatrixDRS(fineOperator, COMPONENT_INDEX, weights_, param)),
scaledMatrixOperator_(Detail::createOperator(fineOperator, *scaledMatrix_, comm)),
smoother_(Detail::constructSmoother<Smoother>(*scaledMatrixOperator_,
smargs, comm)),
levelTransferPolicy_(criterion, comm, param.cpr_pressure_aggregation_),
coarseSolverPolicy_(&param, smargs, criterion),
twoLevelMethod_(*scaledMatrixOperator_,
smoother_,
levelTransferPolicy_,
coarseSolverPolicy_, 0, 1)
{
}
void updatePreconditioner(const Operator& fineOperator,
const SmootherArgs& smargs,
const Communication& comm)
{
*scaledMatrix_ = *Detail::scaleMatrixDRS(fineOperator, COMPONENT_INDEX, weights_, param_);
smoother_ = Detail::constructSmoother<Smoother>(*scaledMatrixOperator_, smargs, comm);
twoLevelMethod_.updatePreconditioner(*scaledMatrixOperator_,
smoother_,
coarseSolverPolicy_);
}
void pre(typename TwoLevelMethod::FineDomainType& x,
typename TwoLevelMethod::FineRangeType& b) override
{
twoLevelMethod_.pre(x,b);
}
void post(typename TwoLevelMethod::FineDomainType& x) override
{
twoLevelMethod_.post(x);
}
void apply(typename TwoLevelMethod::FineDomainType& v,
const typename TwoLevelMethod::FineRangeType& d) override
{
auto scaledD = d;
Detail::scaleVectorDRS(scaledD, COMPONENT_INDEX, param_, weights_);
twoLevelMethod_.apply(v, scaledD);
}
private:
const CPRParameter& param_;
const typename TwoLevelMethod::FineDomainType& weights_;
std::unique_ptr<Matrix> scaledMatrix_;
std::unique_ptr<Operator> scaledMatrixOperator_;
std::shared_ptr<Smoother> smoother_;
LevelTransferPolicy levelTransferPolicy_;
CoarseSolverPolicy coarseSolverPolicy_;
TwoLevelMethod twoLevelMethod_;
};
} // end namespace Opm
#endif

303
opm/simulators/linalg/ISTLSolverEbosCpr.hpp
View File

@ -1,303 +0,0 @@
/*
Copyright 2016 IRIS AS
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_ISTLSOLVERCPR_EBOS_HEADER_INCLUDED
#define OPM_ISTLSOLVERCPR_EBOS_HEADER_INCLUDED
#include <opm/simulators/linalg/ISTLSolverEbos.hpp>
#include <opm/simulators/linalg/BlackoilAmgCpr.hpp>
#include <utility>
#include <memory>
namespace Opm
{
//=====================================================================
// Implementation for ISTL-matrix based operator
//=====================================================================
/// 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
/// 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 ISTLSolverEbosCpr : public ISTLSolverEbos<TypeTag>
{
protected:
// Types and indices from superclass.
using SuperClass = ISTLSolverEbos<TypeTag>;
using Matrix = typename SuperClass::Matrix;
using Vector = typename SuperClass::Vector;
using WellModel = typename SuperClass::WellModel;
using Simulator = typename SuperClass::Simulator;
using SparseMatrixAdapter = typename SuperClass::SparseMatrixAdapter;
enum { pressureEqnIndex = SuperClass::pressureEqnIndex };
enum { pressureVarIndex = SuperClass::pressureVarIndex };
// New properties in this subclass.
using Preconditioner = Dune::Preconditioner<Vector, Vector>;
using MatrixAdapter = Dune::MatrixAdapter<Matrix,Vector, Vector>;
using CouplingMetric = Opm::Amg::Element<pressureEqnIndex,pressureVarIndex>;
using CritBase = Dune::Amg::SymmetricCriterion<Matrix, CouplingMetric>;
using Criterion = Dune::Amg::CoarsenCriterion<CritBase>;
using CprSmootherFine = Opm::ParallelOverlappingILU0<Matrix, Vector, Vector, Dune::Amg::SequentialInformation>;
using CprSmootherCoarse = CprSmootherFine;
using BlackoilAmgType = BlackoilAmgCpr<MatrixAdapter,CprSmootherFine, CprSmootherCoarse, Criterion, Dune::Amg::SequentialInformation,
pressureEqnIndex, pressureVarIndex>;
using OperatorSerial = WellModelMatrixAdapter< Matrix, Vector, Vector, WellModel, false>;
#if HAVE_MPI
using POrComm = Dune::OwnerOverlapCopyCommunication<int,int>;
using ParallelMatrixAdapter = Dune::OverlappingSchwarzOperator<Matrix, Vector, Vector, POrComm >;
using ParallelCprSmootherFine = Opm::ParallelOverlappingILU0<Matrix, Vector, Vector, POrComm >;
using ParallelCprSmootherCoarse = ParallelCprSmootherFine;
using ParallelBlackoilAmgType = BlackoilAmgCpr<ParallelMatrixAdapter, ParallelCprSmootherFine, ParallelCprSmootherCoarse, Criterion,
POrComm, pressureEqnIndex, pressureVarIndex>;
using OperatorParallel = WellModelMatrixAdapter< Matrix, Vector, Vector, WellModel, true>;
using ParallelScalarProduct = Dune::OverlappingSchwarzScalarProduct<Vector, POrComm>;
#else
using POrComm = Dune::Amg::SequentialInformation;
using ParallelBlackoilAmgType = BlackoilAmgType;
using ParallelScalarProduct = Dune::SeqScalarProduct<Vector>;
using ParallelMatrixAdapter = MatrixAdapter;
using OperatorParallel = OperatorSerial;
#endif
public:
static void registerParameters()
{
FlowLinearSolverParameters::registerParameters<TypeTag>();
}
/// Construct a system solver.
/// \param[in] parallelInformation In the case of a parallel run
/// with dune-istl the information about the parallelization.
explicit ISTLSolverEbosCpr(const Simulator& simulator)
: SuperClass(simulator), oldMat()
{}
void prepare(const SparseMatrixAdapter& M, Vector& b)
{
if (oldMat != nullptr)
std::cout << "old was "<<oldMat<<" new is "<<&M.istlMatrix()<<std::endl;
oldMat = &M.istlMatrix();
int newton_iteration = this->simulator_.model().newtonMethod().numIterations();
if (newton_iteration < 1 or not(this->parameters_.cpr_reuse_setup_)) {
SuperClass::matrix_.reset(new Matrix(M.istlMatrix()));
} else {
*SuperClass::matrix_ = M.istlMatrix();
}
SuperClass::rhs_ = &b;
SuperClass::scaleSystem();
const WellModel& wellModel = this->simulator_.problem().wellModel();
#if HAVE_MPI
if( this->isParallel() ) {
//remove ghost rows in local matrix without doing a copy.
this->makeOverlapRowsInvalid(*(this->matrix_));
if (newton_iteration < 1 or not(this->parameters_.cpr_reuse_setup_)) {
//Not sure what actual_mat_for_prec is, so put ebosJacIgnoreOverlap as both variables
//to be certain that correct matrix is used for preconditioning.
if( ! comm_ )
{
opAParallel_.reset(new OperatorParallel(*(this->matrix_), *(this->matrix_), wellModel,
this->parallelInformation_ ));
comm_ = opAParallel_->comm();
assert(comm_->indexSet().size()==0);
const size_t size = opAParallel_->getmat().N();
const ParallelISTLInformation& info =
std::any_cast<const ParallelISTLInformation&>( this->parallelInformation_);
// As we use a dune-istl with block size np the number of components
// per parallel is only one.
info.copyValuesTo(comm_->indexSet(), comm_->remoteIndices(),
size, 1);
}
else
{
opAParallel_.reset(new OperatorParallel(*(this->matrix_), *(this->matrix_), wellModel,
comm_ ));
}
}
constexpr Dune::SolverCategory::Category category=Dune::SolverCategory::overlapping;
auto sp = Dune::createScalarProduct<Vector,POrComm>(*comm_, category);
sp_ = std::move(sp);
using AMGOperator = Dune::OverlappingSchwarzOperator<Matrix, Vector, Vector, POrComm>;
// If clause is always execute as as Linearoperator is WellModelMatrixAdapter< Matrix, Vector, Vector, WellModel, false|true>;
if( ! std::is_same< OperatorParallel, AMGOperator > :: value &&
( newton_iteration < 1 or not(this->parameters_.cpr_reuse_setup_) ) ) {
// create new operator in case linear operator and matrix operator differ
opA_.reset( new AMGOperator( opAParallel_->getmat(), *comm_ ));
}
prepareSolver(*opAParallel_, *comm_);
} else
#endif
{
if (newton_iteration < 1 or not(this->parameters_.cpr_reuse_setup_)) {
opASerial_.reset(new OperatorSerial(*(this->matrix_), *(this->matrix_), wellModel));
}
using POrCommType = Dune::Amg::SequentialInformation;
POrCommType parallelInformation_arg;
typedef OperatorSerial LinearOperator;
constexpr Dune::SolverCategory::Category category=Dune::SolverCategory::sequential;
auto sp = Dune::createScalarProduct<Vector,POrComm>(parallelInformation_arg, category);
sp_ = std::move(sp);
// If clause is always execute as as Linearoperator is WellModelMatrixAdapter< Matrix, Vector, Vector, WellModel, false|true>;
if( ! std::is_same< LinearOperator, MatrixAdapter > :: value &&
( newton_iteration < 1 or not(this->parameters_.cpr_reuse_setup_) ) ) {
// create new operator in case linear operator and matrix operator differ
opA_.reset( new MatrixAdapter( opASerial_->getmat()));//, parallelInformation_arg ) );
}
prepareSolver(*opASerial_, parallelInformation_arg);
}
}
template<typename Operator, typename Comm>
void prepareSolver(Operator& wellOpA, Comm& comm)
{
Vector& istlb = *(this->rhs_);
comm.copyOwnerToAll(istlb, istlb);
const double relax = this->parameters_.ilu_relaxation_;
const MILU_VARIANT ilu_milu = this->parameters_.ilu_milu_;
// TODO: revise choice of parameters
// int coarsenTarget = 4000;
int coarsenTarget = 1200;
Criterion criterion(15, coarsenTarget);
criterion.setDebugLevel( this->parameters_.cpr_solver_verbose_ ); // no debug information, 1 for printing hierarchy information
criterion.setDefaultValuesIsotropic(2);
criterion.setNoPostSmoothSteps( 1 );
criterion.setNoPreSmoothSteps( 1 );
//new guesses by hmbn
//criterion.setAlpha(0.01); // criterion for connection strong 1/3 is default
//criterion.setMaxLevel(2); //
//criterion.setGamma(1); // //1 V cycle 2 WW
// Since DUNE 2.2 we also need to pass the smoother args instead of steps directly
using AmgType = typename std::conditional<std::is_same<Comm, Dune::Amg::SequentialInformation>::value,
BlackoilAmgType, ParallelBlackoilAmgType>::type;
using SpType = typename std::conditional<std::is_same<Comm, Dune::Amg::SequentialInformation>::value,
Dune::SeqScalarProduct<Vector>,
ParallelScalarProduct >::type;
using OperatorType = typename std::conditional<std::is_same<Comm, Dune::Amg::SequentialInformation>::value,
MatrixAdapter, ParallelMatrixAdapter>::type;
typedef typename AmgType::Smoother Smoother;
typedef typename Dune::Amg::SmootherTraits<Smoother>::Arguments SmootherArgs;
SmootherArgs smootherArgs;
smootherArgs.iterations = 1;
smootherArgs.relaxationFactor = relax;
const Opm::CPRParameter& params(this->parameters_); // strange conversion
ISTLUtility::setILUParameters(smootherArgs, ilu_milu);
auto& opARef = reinterpret_cast<OperatorType&>(*opA_);
int newton_iteration = this->simulator_.model().newtonMethod().numIterations();
bool update_preconditioner = false;
if (this->parameters_.cpr_reuse_setup_ < 1) {
update_preconditioner = true;
}
if (this->parameters_.cpr_reuse_setup_ < 2) {
if (newton_iteration < 1) {
update_preconditioner = true;
}
}
if (this->parameters_.cpr_reuse_setup_ < 3) {
if (this->iterations() > 10) {
update_preconditioner = true;
}
}
if ( update_preconditioner or (amg_== 0) ) {
amg_.reset( new AmgType( params, this->weights_, opARef, criterion, smootherArgs, comm ) );
} else {
if (this->parameters_.cpr_solver_verbose_) {
std::cout << " Only update amg solver " << std::endl;
}
reinterpret_cast<AmgType*>(amg_.get())->updatePreconditioner(opARef, smootherArgs, comm);
}
// Solve.
//SuperClass::solve(linearOperator, x, istlb, *sp, *amg, result);
//references seems to do something els than refering
int verbosity_linsolve = 0;
if (comm.communicator().rank() == 0) {
verbosity_linsolve = this->parameters_.linear_solver_verbosity_;
}
linsolve_.reset(new Dune::BiCGSTABSolver<Vector>(wellOpA, reinterpret_cast<SpType&>(*sp_), reinterpret_cast<AmgType&>(*amg_),
this->parameters_.linear_solver_reduction_,
this->parameters_.linear_solver_maxiter_,
verbosity_linsolve));
}
bool solve(Vector& x)
{
// Solve system.
Dune::InverseOperatorResult result;
Vector& istlb = *(this->rhs_);
linsolve_->apply(x, istlb, result);
SuperClass::checkConvergence(result);
if (this->parameters_.scale_linear_system_) {
this->scaleSolution(x);
}
return this->converged_;
}
protected:
///! \brief The dune-istl operator (either serial or parallel
std::unique_ptr< Dune::LinearOperator<Vector, Vector> > opA_;
///! \brief Serial well matrix adapter
std::unique_ptr< OperatorSerial > opASerial_;
///! \brief Parallel well matrix adapter
std::unique_ptr< OperatorParallel > opAParallel_;
///! \brief The preconditoner to use (either serial or parallel CPR with AMG)
std::unique_ptr< Preconditioner > amg_;
using SPPointer = std::shared_ptr< Dune::ScalarProduct<Vector> >;
SPPointer sp_;
std::shared_ptr< Dune::BiCGSTABSolver<Vector> > linsolve_;
const void* oldMat;
std::shared_ptr<POrComm> comm_;
}; // end ISTLSolver
} // namespace Opm
#endif