mirror of
https://github.com/OPM/opm-simulators.git
synced 2026-10-12 08:54:53 -05:00
Clean up formatting and get most types from superclass.
This commit is contained in:
Atgeirr Flø Rasmussen authored and 
Markus Blatt committed
1 parent
be86b90bac
commit
7664d27849
1 file changed
+164
-239
+164
-239
@@ -48,53 +48,37 @@ namespace Opm
|
||||
template <class TypeTag>
|
||||
class ISTLSolverEbosCpr : public ISTLSolverEbos<TypeTag>
|
||||
{
|
||||
typedef typename GET_PROP_TYPE(TypeTag, GridView) GridView;
|
||||
typedef typename GET_PROP_TYPE(TypeTag, Scalar) Scalar;
|
||||
typedef typename GET_PROP_TYPE(TypeTag, SparseMatrixAdapter) SparseMatrixAdapter;
|
||||
typedef typename GET_PROP_TYPE(TypeTag, GlobalEqVector) Vector;
|
||||
typedef typename GET_PROP_TYPE(TypeTag, Indices) Indices;
|
||||
typedef typename GET_PROP_TYPE(TypeTag, EclWellModel) WellModel;
|
||||
typedef typename GET_PROP_TYPE(TypeTag, Simulator) Simulator;
|
||||
typedef typename SparseMatrixAdapter::IstlMatrix Matrix;
|
||||
typedef typename GET_PROP_TYPE(TypeTag, CprSmootherFine) CprSmootherFine;
|
||||
typedef typename GET_PROP_TYPE(TypeTag, CprSmootherCoarse) CprSmootherCoarse;
|
||||
//typedef typename SparseMatrixAdapter::MatrixBlock MatrixBlockType;
|
||||
//typedef typename Vector::block_type BlockVector;
|
||||
//typedef typename GET_PROP_TYPE(TypeTag, Evaluation) Evaluation;
|
||||
//typedef typename GET_PROP_TYPE(TypeTag, ThreadManager) ThreadManager;
|
||||
//typedef typename GridView::template Codim<0>::Entity Element;
|
||||
//typedef typename GET_PROP_TYPE(TypeTag, ElementContext) ElementContext;
|
||||
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 };
|
||||
|
||||
enum { pressureEqnIndex = BlackOilDefaultIndexTraits::waterCompIdx };
|
||||
enum { pressureVarIndex = Indices::pressureSwitchIdx };
|
||||
// New properties in this subclass.
|
||||
using CprSmootherFine = typename GET_PROP_TYPE(TypeTag, CprSmootherFine);
|
||||
using CprSmootherCoarse = typename GET_PROP_TYPE(TypeTag, CprSmootherCoarse);
|
||||
|
||||
using OperatorSerial = WellModelMatrixAdapter< Matrix, Vector, Vector, WellModel, false>;
|
||||
using POrComm = Dune::Amg::SequentialInformation;
|
||||
using MatrixAdapter = Dune::MatrixAdapter<Matrix,Vector, Vector>;
|
||||
|
||||
#if DUNE_VERSION_NEWER(DUNE_ISTL, 2, 6)
|
||||
|
||||
static const int numEq = Indices::numEq;
|
||||
typedef WellModelMatrixAdapter< Matrix, Vector, Vector, WellModel, false> OperatorSerial;
|
||||
typedef ISTLSolverEbos<TypeTag> SuperClass;
|
||||
typedef Dune::Amg::SequentialInformation POrComm;
|
||||
//typedef ISTLUtility::CPRSelector< Matrix, Vector, Vector, POrComm> CPRSelectorType;
|
||||
typedef Dune::MatrixAdapter<Matrix,Vector, Vector> MatrixAdapter;
|
||||
|
||||
#if DUNE_VERSION_NEWER(DUNE_ISTL, 2, 6)
|
||||
|
||||
#else
|
||||
static constexpr int category = Dune::SolverCategory::sequential;
|
||||
typedef Dune::ScalarProductChooser<Vector, POrComm, category> ScalarProductChooser;
|
||||
#endif
|
||||
//Operator MatrixOperator = Dune::MatrixAdapter<Matrix,Vector,Vector>
|
||||
//typedef Opm::ParallelOverlappingILU0<Matrix,Vector,Vector, POrComm> Smoother;
|
||||
typedef CprSmootherFine Smoother;
|
||||
//ParallelInformation = Dune::Amg::SequentialInformation
|
||||
//typedef Dune::Amg::AMG<MatrixAdapter,Vector,Smoother,POrComm> DuneAmg;
|
||||
using CouplingMetric = Opm::Amg::Element<pressureEqnIndex,pressureVarIndex>;
|
||||
using CritBase = Dune::Amg::SymmetricCriterion<Matrix, CouplingMetric>;
|
||||
using Criterion = Dune::Amg::CoarsenCriterion<CritBase>;
|
||||
typedef BlackoilAmgCpr<MatrixAdapter,CprSmootherFine, CprSmootherCoarse, Criterion, POrComm, pressureEqnIndex, pressureVarIndex> BlackoilAmgType;
|
||||
|
||||
using CouplingMetric = Opm::Amg::Element<pressureEqnIndex,pressureVarIndex>;
|
||||
using CritBase = Dune::Amg::SymmetricCriterion<Matrix, CouplingMetric>;
|
||||
using Criterion = Dune::Amg::CoarsenCriterion<CritBase>;
|
||||
using BlackoilAmgType = BlackoilAmgCpr<MatrixAdapter,CprSmootherFine, CprSmootherCoarse, Criterion, POrComm, pressureEqnIndex, pressureVarIndex>;
|
||||
|
||||
public:
|
||||
typedef Dune::AssembledLinearOperator< Matrix, Vector, Vector > AssembledLinearOperatorType;
|
||||
|
||||
static void registerParameters()
|
||||
{
|
||||
FlowLinearSolverParameters::registerParameters<TypeTag>();
|
||||
@@ -103,236 +87,177 @@ namespace Opm
|
||||
/// Construct a system solver.
|
||||
/// \param[in] parallelInformation In the case of a parallel run
|
||||
/// with dune-istl the information about the parallelization.
|
||||
ISTLSolverEbosCpr(const Simulator& simulator)
|
||||
: SuperClass(simulator)
|
||||
{
|
||||
}
|
||||
|
||||
// nothing to clean here
|
||||
void eraseMatrix() {
|
||||
this->matrix_for_preconditioner_.reset();
|
||||
explicit ISTLSolverEbosCpr(const Simulator& simulator)
|
||||
: SuperClass(simulator)
|
||||
{
|
||||
}
|
||||
|
||||
void prepare(const SparseMatrixAdapter& M, Vector& b){
|
||||
int newton_iteration = this->simulator_.model().newtonMethod().numIterations();
|
||||
// double dt = this->simulator_.timeStepSize();
|
||||
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();
|
||||
void prepare(const SparseMatrixAdapter& M, Vector& b)
|
||||
{
|
||||
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() )
|
||||
{
|
||||
// parallel implemantation si as before
|
||||
// typedef WellModelMatrixAdapter< Matrix, Vector, Vector, WellModel, true ,TypeTag> Operator;
|
||||
if( this->isParallel() ) {
|
||||
// parallel implemantation si as before
|
||||
// typedef WellModelMatrixAdapter< Matrix, Vector, Vector, WellModel, true ,TypeTag> Operator;
|
||||
|
||||
// auto ebosJacIgnoreOverlap = Matrix(*(this->matrix_));
|
||||
// //remove ghost rows in local matrix
|
||||
// this->makeOverlapRowsInvalid(ebosJacIgnoreOverlap);
|
||||
// auto ebosJacIgnoreOverlap = Matrix(*(this->matrix_));
|
||||
// //remove ghost rows in local matrix
|
||||
// this->makeOverlapRowsInvalid(ebosJacIgnoreOverlap);
|
||||
|
||||
// //Not sure what actual_mat_for_prec is, so put ebosJacIgnoreOverlap as both variables
|
||||
// //to be certain that correct matrix is used for preconditioning.
|
||||
// Operator opA(ebosJacIgnoreOverlap, ebosJacIgnoreOverlap, wellModel,
|
||||
// this->parallelInformation_ );
|
||||
// assert( opA.comm() );
|
||||
// //SuperClass::solve( opA, x, *(this->rhs_), *(opA.comm()) );
|
||||
// typedef Dune::OwnerOverlapCopyCommunication<int,int>& comm = *(opA.comm());
|
||||
// const size_t size = opA.getmat().N();
|
||||
// //Not sure what actual_mat_for_prec is, so put ebosJacIgnoreOverlap as both variables
|
||||
// //to be certain that correct matrix is used for preconditioning.
|
||||
// Operator opA(ebosJacIgnoreOverlap, ebosJacIgnoreOverlap, wellModel,
|
||||
// this->parallelInformation_ );
|
||||
// assert( opA.comm() );
|
||||
// //SuperClass::solve( opA, x, *(this->rhs_), *(opA.comm()) );
|
||||
// typedef Dune::OwnerOverlapCopyCommunication<int,int>& comm = *(opA.comm());
|
||||
// const size_t size = opA.getmat().N();
|
||||
|
||||
// const ParallelISTLInformation& info =
|
||||
// boost::any_cast<const ParallelISTLInformation&>( this->parallelInformation_);
|
||||
// const ParallelISTLInformation& info =
|
||||
// boost::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);
|
||||
// // Construct operator, scalar product and vectors needed.
|
||||
// Dune::InverseOperatorResult result;
|
||||
// SuperClass::constructPreconditionerAndSolve<Dune::SolverCategory::overlapping>(opA, x, *(this->rhs_), comm, result);
|
||||
// SuperClass::checkConvergence(result);
|
||||
// // 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);
|
||||
// // Construct operator, scalar product and vectors needed.
|
||||
// Dune::InverseOperatorResult result;
|
||||
// SuperClass::constructPreconditionerAndSolve<Dune::SolverCategory::overlapping>(opA, x, *(this->rhs_), comm, result);
|
||||
// SuperClass::checkConvergence(result);
|
||||
|
||||
}
|
||||
else
|
||||
} else
|
||||
#endif
|
||||
{
|
||||
const WellModel& wellModel = this->simulator_.problem().wellModel();
|
||||
//typedef WellModelMatrixAdapter< Matrix, Vector, Vector, WellModel, false ,TypeTag> OperatorSerial;
|
||||
opASerial_.reset(new OperatorSerial(*(this->matrix_), *(this->matrix_), wellModel));
|
||||
|
||||
//Dune::Amg::SequentialInformation info;
|
||||
typedef Dune::Amg::SequentialInformation POrComm;
|
||||
POrComm parallelInformation_arg;
|
||||
typedef OperatorSerial LinearOperator;
|
||||
|
||||
//SuperClass::constructPreconditionerAndSolve(opA, x, *(this->rhs_), info, result);
|
||||
|
||||
|
||||
{
|
||||
const WellModel& wellModel = this->simulator_.problem().wellModel();
|
||||
opASerial_.reset(new OperatorSerial(*(this->matrix_), *(this->matrix_), wellModel));
|
||||
typedef Dune::Amg::SequentialInformation POrComm;
|
||||
POrComm parallelInformation_arg;
|
||||
typedef OperatorSerial LinearOperator;
|
||||
|
||||
#if DUNE_VERSION_NEWER(DUNE_ISTL, 2, 6)
|
||||
constexpr Dune::SolverCategory::Category category=Dune::SolverCategory::sequential;
|
||||
auto sp = Dune::createScalarProduct<Vector,POrComm>(parallelInformation_arg, category);
|
||||
sp_ = std::move(sp);
|
||||
constexpr Dune::SolverCategory::Category category=Dune::SolverCategory::sequential;
|
||||
auto sp = Dune::createScalarProduct<Vector,POrComm>(parallelInformation_arg, category);
|
||||
sp_ = std::move(sp);
|
||||
#else
|
||||
constexpr int category = Dune::SolverCategory::sequential;
|
||||
typedef Dune::ScalarProductChooser<Vector, POrComm, category> ScalarProductChooser;
|
||||
typedef std::unique_ptr<typename ScalarProductChooser::ScalarProduct> SPPointer;
|
||||
SPPointer sp(ScalarProductChooser::construct(parallelInformation_arg));
|
||||
sp_ = std::move(sp);
|
||||
constexpr int category = Dune::SolverCategory::sequential;
|
||||
typedef Dune::ScalarProductChooser<Vector, POrComm, category> ScalarProductChooser;
|
||||
typedef std::unique_ptr<typename ScalarProductChooser::ScalarProduct> SPPointer;
|
||||
SPPointer sp(ScalarProductChooser::construct(parallelInformation_arg));
|
||||
sp_ = std::move(sp);
|
||||
#endif
|
||||
Vector& istlb = *(this->rhs_);
|
||||
parallelInformation_arg.copyOwnerToAll(istlb, istlb);
|
||||
|
||||
|
||||
Vector& istlb = *(this->rhs_);
|
||||
parallelInformation_arg.copyOwnerToAll(istlb, istlb);
|
||||
|
||||
if( ! std::is_same< LinearOperator, MatrixAdapter > :: value )
|
||||
{
|
||||
// create new operator in case linear operator and matrix operator differ
|
||||
opA_.reset( new MatrixAdapter( opASerial_->getmat()));//, parallelInformation_arg ) );
|
||||
if( ! std::is_same< LinearOperator, MatrixAdapter > :: value ) {
|
||||
// create new operator in case linear operator and matrix operator differ
|
||||
opA_.reset( new MatrixAdapter( opASerial_->getmat()));//, parallelInformation_arg ) );
|
||||
}
|
||||
|
||||
const double relax = this->parameters_.ilu_relaxation_;
|
||||
const MILU_VARIANT ilu_milu = this->parameters_.ilu_milu_;
|
||||
using Matrix = typename MatrixAdapter::matrix_type;
|
||||
//using CouplingMetric = Dune::Amg::Diagonal<pressureIndex>;
|
||||
//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;
|
||||
|
||||
//std::unique_ptr< AMG > amg;
|
||||
// Construct preconditioner.
|
||||
//Criterion crit(15, 2000);
|
||||
//SuperClass::constructAMGPrecond< Criterion >( linearOperator, parallelInformation_arg, amg, opA, relax, ilu_milu );
|
||||
// ISTLUtility::template createAMGPreconditionerPointer<Criterion>( *opA_,
|
||||
// relax,
|
||||
// parallelInformation_arg,
|
||||
// amg_,
|
||||
// this->parameters_,
|
||||
// this->weights_ );
|
||||
//using AMG = BlackoilAmg<Op,S,Criterion,P, PressureIndex>;
|
||||
POrComm& comm = parallelInformation_arg;
|
||||
const int verbosity = ( this->parameters_.cpr_solver_verbose_ &&
|
||||
comm.communicator().rank()==0 ) ? 1 : 0;
|
||||
const double relax = this->parameters_.ilu_relaxation_;
|
||||
const MILU_VARIANT ilu_milu = this->parameters_.ilu_milu_;
|
||||
using Matrix = typename MatrixAdapter::matrix_type;
|
||||
POrComm& comm = parallelInformation_arg;
|
||||
const int verbosity = ( this->parameters_.cpr_solver_verbose_ &&
|
||||
comm.communicator().rank()==0 ) ? 1 : 0;
|
||||
|
||||
// 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
|
||||
// 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
|
||||
typedef typename BlackoilAmgType::Smoother Smoother;
|
||||
typedef typename BlackoilAmgType::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);
|
||||
//ISTLUtility::setILUParameters(smootherArgs, params);
|
||||
//smootherArgs.setN(params.cpr_ilu_n_); smootherArgs.setMilu(params.cpr_ilu_milu_);
|
||||
// Since DUNE 2.2 we also need to pass the smoother args instead of steps directly
|
||||
typedef typename BlackoilAmgType::Smoother Smoother;
|
||||
typedef typename BlackoilAmgType::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);
|
||||
|
||||
MatrixAdapter& opARef = *opA_;
|
||||
int newton_iteration = this->simulator_.model().newtonMethod().numIterations();
|
||||
double dt = this->simulator_.timeStepSize();
|
||||
bool update_preconditioner = false;
|
||||
MatrixAdapter& opARef = *opA_;
|
||||
int newton_iteration = this->simulator_.model().newtonMethod().numIterations();
|
||||
double dt = this->simulator_.timeStepSize();
|
||||
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 (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 BlackoilAmgType( params, this->weights_, opARef, criterion, smootherArgs, comm ) );
|
||||
}else{
|
||||
if(this->parameters_.cpr_solver_verbose_){
|
||||
std::cout << " Only update amg solver " << std::endl;
|
||||
}
|
||||
amg_->updatePreconditioner(opARef, smootherArgs, comm);
|
||||
}
|
||||
// Solve.
|
||||
//SuperClass::solve(linearOperator, x, istlb, *sp, *amg, result);
|
||||
//references seems to do something els than refering
|
||||
if ( update_preconditioner or (amg_== 0) ) {
|
||||
amg_.reset( new BlackoilAmgType( params, this->weights_, opARef, criterion, smootherArgs, comm ) );
|
||||
} else {
|
||||
if (this->parameters_.cpr_solver_verbose_) {
|
||||
std::cout << " Only update amg solver " << std::endl;
|
||||
}
|
||||
amg_->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_;
|
||||
}
|
||||
int verbosity_linsolve = 0;
|
||||
if (comm.communicator().rank() == 0) {
|
||||
verbosity_linsolve = this->parameters_.linear_solver_verbosity_;
|
||||
}
|
||||
|
||||
LinearOperator& opASerialRef = *opASerial_;
|
||||
linsolve_.reset(new Dune::BiCGSTABSolver<Vector>(opASerialRef, *sp_, *amg_,
|
||||
this->parameters_.linear_solver_reduction_,
|
||||
this->parameters_.linear_solver_maxiter_,
|
||||
verbosity_linsolve));
|
||||
|
||||
|
||||
}
|
||||
LinearOperator& opASerialRef = *opASerial_;
|
||||
linsolve_.reset(new Dune::BiCGSTABSolver<Vector>(opASerialRef, *sp_, *amg_,
|
||||
this->parameters_.linear_solver_reduction_,
|
||||
this->parameters_.linear_solver_maxiter_,
|
||||
verbosity_linsolve));
|
||||
}
|
||||
}
|
||||
|
||||
bool solve(Vector& x) {
|
||||
//SuperClass::solve(x);
|
||||
if( this->isParallel() ){
|
||||
// for now only call the superclass
|
||||
bool converged = SuperClass::solve(x);
|
||||
return converged;
|
||||
}else{
|
||||
// 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);
|
||||
}
|
||||
}
|
||||
bool solve(Vector& x)
|
||||
{
|
||||
if (this->isParallel()) {
|
||||
// for now only call the superclass
|
||||
bool converged = SuperClass::solve(x);
|
||||
return converged;
|
||||
} else {
|
||||
// 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_;
|
||||
}
|
||||
|
||||
|
||||
/// Solve the system of linear equations Ax = b, with A being the
|
||||
/// combined derivative matrix of the residual and b
|
||||
/// being the residual itself.
|
||||
/// \param[in] residual residual object containing A and b.
|
||||
/// \return the solution x
|
||||
|
||||
/// \copydoc NewtonIterationBlackoilInterface::iterations
|
||||
int iterations () const { return this->iterations_; }
|
||||
|
||||
/// \copydoc NewtonIterationBlackoilInterface::parallelInformation
|
||||
const boost::any& parallelInformation() const { return this->parallelInformation_; }
|
||||
|
||||
protected:
|
||||
|
||||
//using Matrix = typename MatrixAdapter::matrix_type;
|
||||
//using CouplingMetric = Dune::Amg::Diagonal<pressureIndex>;
|
||||
//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;
|
||||
//Operator MatrixOperator = Dune::MatrixAdapter<Matrix,Vector,Vector>
|
||||
//Smoother ParallelOverLappingILU0<Matrix,Vector,Vector>
|
||||
//ParallelInformation = Dune::Amg::SequentialInformation
|
||||
|
||||
#if DUNE_VERSION_NEWER(DUNE_ISTL, 2, 6)
|
||||
typedef std::shared_ptr< Dune::ScalarProduct<Vector> > SPPointer;
|
||||
#else
|
||||
|
||||
Reference in new issue
Block a user