mirror of
https://github.com/OPM/opm-simulators.git
synced 2025-02-25 18:55:30 -06:00
1c2d3fbcc7
That the local ids were consecutive and starting from 0 was just a coincidence and they should never be used to access linear systems or vectors. This commit fixes this by using the correct mappers instead. Note the we removed some computations from the constructor of ISTLSolverEbosCpr as it inherits from ISTLSolverEbos and the operations already happnen in constructor of the base class.
304 lines
15 KiB
C++
304 lines
15 KiB
C++
/*
|
|
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
|