OilfieldToolsNet/opm-simulators
4
0
Fork 0
mirror of https://github.com/OPM/opm-simulators.git synced 2024-12-23 07:53:29 -06:00
Code Issues Packages Projects Releases Wiki Activity

Adapt to linear solver api from ewoms

Browse Source
This commit is contained in:
Tor Harald Sandve 2018-11-12 11:03:54 +01:00
parent 578923f93a
commit 2d2d0c4433
8 changed files with 241 additions and 239 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

1
opm/autodiff/BlackoilDetails.hpp
View File

@ -25,6 +25,7 @@
#define OPM_BLACKOILDETAILS_HEADER_INCLUDED
#include <opm/core/linalg/ParallelIstlInformation.hpp>
#include <opm/core/props/BlackoilPhases.hpp>
namespace Opm {
namespace detail {

200
opm/autodiff/BlackoilModelEbos.hpp
View File

@ -86,6 +86,8 @@ SET_BOOL_PROP(EclFlowProblem, EnableTemperature, true);
SET_BOOL_PROP(EclFlowProblem, EnableEnergy, false);
SET_TYPE_PROP(EclFlowProblem, EclWellModel, Opm::BlackoilWellModel<TypeTag>);
SET_TAG_PROP(EclFlowProblem, LinearSolverSplice, FlowIstlSolver);
END_PROPERTIES
@ -149,11 +151,9 @@ namespace Opm {
BlackoilModelEbos(Simulator& ebosSimulator,
const ModelParameters& param,
BlackoilWellModel<TypeTag>& well_model,
const ISTLSolverType& linsolver,
const bool terminal_output)
: ebosSimulator_(ebosSimulator)
, grid_(ebosSimulator_.vanguard().grid())
, istlSolver_( dynamic_cast< const ISTLSolverType* > (&linsolver) )
, phaseUsage_(phaseUsageFromDeck(eclState()))
, has_disgas_(FluidSystem::enableDissolvedGas())
, has_vapoil_(FluidSystem::enableVaporizedOil())
@ -169,12 +169,6 @@ namespace Opm {
{
// compute global sum of number of cells
global_nc_ = detail::countGlobalCells(grid_);
//find rows of matrix corresponding to overlap
detail::findOverlapRowsAndColumns(grid_,overlapRowAndColumns_);
if (!istlSolver_)
{
OPM_THROW(std::logic_error,"solver down cast to ISTLSolver failed");
}
convergence_reports_.reserve(300); // Often insufficient, but avoids frequent moves.
}
@ -375,16 +369,6 @@ namespace Opm {
ebosSimulator_.model().linearizer().linearize();
ebosSimulator_.problem().endIteration();
auto& ebosJac = ebosSimulator_.model().linearizer().jacobian();
if (param_.matrix_add_well_contributions_) {
wellModel().addWellContributions(ebosJac.istlMatrix());
}
if ( param_.preconditioner_add_well_contributions_ &&
! param_.matrix_add_well_contributions_ ) {
matrix_for_preconditioner_ .reset(new Mat(ebosJac.istlMatrix()));
wellModel().addWellContributions(*matrix_for_preconditioner_);
}
return wellModel().lastReport();
}
@ -469,179 +453,39 @@ namespace Opm {
/// Number of linear iterations used in last call to solveJacobianSystem().
int linearIterationsLastSolve() const
{
return istlSolver().iterations();
}
/// Zero out off-diagonal blocks on rows corresponding to overlap cells
/// Diagonal blocks on ovelap rows are set to diag(1e100).
void makeOverlapRowsInvalid(Mat& ebosJacIgnoreOverlap) const
{
//value to set on diagonal
MatrixBlockType diag_block(0.0);
for (int eq = 0; eq < numEq; ++eq)
diag_block[eq][eq] = 1.0e100;
//loop over precalculated overlap rows and columns
for (auto row = overlapRowAndColumns_.begin(); row != overlapRowAndColumns_.end(); row++ )
{
int lcell = row->first;
//diagonal block set to large value diagonal
ebosJacIgnoreOverlap[lcell][lcell] = diag_block;
//loop over off diagonal blocks in overlap row
for (auto col = row->second.begin(); col != row->second.end(); ++col)
{
int ncell = *col;
//zero out block
ebosJacIgnoreOverlap[lcell][ncell] = 0.0;
}
}
return ebosSimulator_.model().newtonMethod().linearSolver().iterations ();
}
/// Solve the Jacobian system Jx = r where J is the Jacobian and
/// r is the residual.
void solveJacobianSystem(BVector& x) const
void solveJacobianSystem(BVector& x)
{
auto& ebosJac = ebosSimulator_.model().linearizer().jacobian();
auto& ebosResid = ebosSimulator_.model().linearizer().residual();
// J = [A, B; C, D], where A is the reservoir equations, B and C the interaction of well
// with the reservoir and D is the wells itself.
// The full system is reduced to a number of cells X number of cells system via Schur complement
// A -= B^T D^-1 C
// Instead of modifying A, the Ax operator is modified. i.e Ax -= B^T D^-1 C x in the WellModelMatrixAdapter.
// If matrix_add_well_contribution is false, the Ax operator is modified. i.e Ax -= B^T D^-1 C x in the WellModelMatrixAdapter
// instead of A.
// The residual is modified similarly.
// r = [r, r_well], where r is the residual and r_well the well residual.
// r -= B^T * D^-1 r_well
auto& ebosJac = ebosSimulator_.model().linearizer().jacobian();
auto& ebosResid = ebosSimulator_.model().linearizer().residual();
wellModel().apply(ebosResid);
if (param_.matrix_add_well_contributions_) {
wellModel().addWellContributions(ebosJac.istlMatrix());
}
// set initial guess
x = 0.0;
const Mat& actual_mat_for_prec = matrix_for_preconditioner_ ? *matrix_for_preconditioner_.get() : ebosJac.istlMatrix();
// Solve system.
if( isParallel() )
{
typedef WellModelMatrixAdapter< Mat, BVector, BVector, BlackoilWellModel<TypeTag>, true > Operator;
auto& ebosSolver = ebosSimulator_.model().newtonMethod().linearSolver();
ebosSolver.prepare(ebosJac, ebosResid);
ebosSolver.solve(x);
}
auto ebosJacIgnoreOverlap = Mat(ebosJac.istlMatrix());
//remove ghost rows in local matrix
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(),
istlSolver().parallelInformation() );
assert( opA.comm() );
istlSolver().solve( opA, x, ebosResid, *(opA.comm()) );
}
else
{
typedef WellModelMatrixAdapter< Mat, BVector, BVector, BlackoilWellModel<TypeTag>, false > Operator;
Operator opA(ebosJac.istlMatrix(), actual_mat_for_prec, wellModel());
istlSolver().solve( opA, x, ebosResid );
}
}
//=====================================================================
// Implementation for ISTL-matrix based operator
//=====================================================================
/*!
\brief Adapter to turn a matrix into a linear operator.
Adapts a matrix to the assembled linear operator interface
*/
template<class M, class X, class Y, class WellModel, bool overlapping >
class WellModelMatrixAdapter : public Dune::AssembledLinearOperator<M,X,Y>
{
typedef Dune::AssembledLinearOperator<M,X,Y> BaseType;
public:
typedef M matrix_type;
typedef X domain_type;
typedef Y range_type;
typedef typename X::field_type field_type;
#if HAVE_MPI
typedef Dune::OwnerOverlapCopyCommunication<int,int> communication_type;
#else
typedef Dune::CollectiveCommunication< Grid > communication_type;
#endif
#if DUNE_VERSION_NEWER(DUNE_ISTL, 2, 6)
Dune::SolverCategory::Category category() const override
{
return overlapping ?
Dune::SolverCategory::overlapping : Dune::SolverCategory::sequential;
}
#else
enum {
//! \brief The solver category.
category = overlapping ?
Dune::SolverCategory::overlapping :
Dune::SolverCategory::sequential
};
#endif
//! constructor: just store a reference to a matrix
WellModelMatrixAdapter (const M& A,
const M& A_for_precond,
const WellModel& wellMod,
const boost::any& parallelInformation = boost::any() )
: A_( A ), A_for_precond_(A_for_precond), wellMod_( wellMod ), comm_()
{
#if HAVE_MPI
if( parallelInformation.type() == typeid(ParallelISTLInformation) )
{
const ParallelISTLInformation& info =
boost::any_cast<const ParallelISTLInformation&>( parallelInformation);
comm_.reset( new communication_type( info.communicator() ) );
}
#endif
}
virtual void apply( const X& x, Y& y ) const
{
A_.mv( x, y );
// add well model modification to y
wellMod_.apply(x, y );
#if HAVE_MPI
if( comm_ )
comm_->project( y );
#endif
}
// y += \alpha * A * x
virtual void applyscaleadd (field_type alpha, const X& x, Y& y) const
{
A_.usmv(alpha,x,y);
// add scaled well model modification to y
wellMod_.applyScaleAdd( alpha, x, y );
#if HAVE_MPI
if( comm_ )
comm_->project( y );
#endif
}
virtual const matrix_type& getmat() const { return A_for_precond_; }
communication_type* comm()
{
return comm_.operator->();
}
protected:
const matrix_type& A_ ;
const matrix_type& A_for_precond_ ;
const WellModel& wellMod_;
std::unique_ptr< communication_type > comm_;
};
/// Apply an update to the primary variables.
void updateSolution(const BVector& dx)
@ -945,6 +789,15 @@ namespace Opm {
std::vector<Scalar> B_avg(numEq, 0.0);
auto report = getReservoirConvergence(timer.currentStepLength(), iteration, B_avg, residual_norms);
report += wellModel().getWellConvergence(B_avg);
// Throw if any NaN or too large residual found.
ConvergenceReport::Severity severity = report.severityOfWorstFailure();
if (severity == ConvergenceReport::Severity::NotANumber) {
OPM_THROW(Opm::NumericalIssue, "NaN residual found!");
} else if (severity == ConvergenceReport::Severity::TooLarge) {
OPM_THROW(Opm::NumericalIssue, "Too large residual found!");
}
return report;
}
@ -1030,9 +883,6 @@ namespace Opm {
double current_relaxation_;
BVector dx_old_;
std::unique_ptr<Mat> matrix_for_preconditioner_;
std::vector<std::pair<int,std::vector<int>>> overlapRowAndColumns_;
std::vector<StepReport> convergence_reports_;
public:
/// return the StandardWells object

7
opm/autodiff/BlackoilModelParametersEbos.hpp
View File

@ -48,7 +48,6 @@ NEW_PROP_TAG(SolveWelleqInitially);
NEW_PROP_TAG(UpdateEquationsScaling);
NEW_PROP_TAG(UseUpdateStabilization);
NEW_PROP_TAG(MatrixAddWellContributions);
NEW_PROP_TAG(PreconditionerAddWellContributions);
// parameters for multisegment wells
NEW_PROP_TAG(TolerancePressureMsWells);
@ -72,7 +71,6 @@ SET_BOOL_PROP(FlowModelParameters, SolveWelleqInitially, true);
SET_BOOL_PROP(FlowModelParameters, UpdateEquationsScaling, false);
SET_BOOL_PROP(FlowModelParameters, UseUpdateStabilization, true);
SET_BOOL_PROP(FlowModelParameters, MatrixAddWellContributions, false);
SET_BOOL_PROP(FlowModelParameters, PreconditionerAddWellContributions, false);
SET_SCALAR_PROP(FlowModelParameters, TolerancePressureMsWells, 0.01 *1e5);
SET_SCALAR_PROP(FlowModelParameters, MaxPressureChangeMsWells, 2.0 *1e5);
SET_BOOL_PROP(FlowModelParameters, UseInnerIterationsMsWells, true);
@ -155,9 +153,6 @@ namespace Opm
// Whether to add influences of wells between cells to the matrix and preconditioner matrix
bool matrix_add_well_contributions_;
// Whether to add influences of wells between cells to the preconditioner matrix only
bool preconditioner_add_well_contributions_;
/// Construct from user parameters or defaults.
BlackoilModelParametersEbos()
{
@ -181,7 +176,6 @@ namespace Opm
update_equations_scaling_ = EWOMS_GET_PARAM(TypeTag, bool, UpdateEquationsScaling);
use_update_stabilization_ = EWOMS_GET_PARAM(TypeTag, bool, UseUpdateStabilization);
matrix_add_well_contributions_ = EWOMS_GET_PARAM(TypeTag, bool, MatrixAddWellContributions);
preconditioner_add_well_contributions_ = EWOMS_GET_PARAM(TypeTag, bool, PreconditionerAddWellContributions);
deck_file_name_ = EWOMS_GET_PARAM(TypeTag, std::string, EclDeckFileName);
}
@ -208,7 +202,6 @@ namespace Opm
EWOMS_REGISTER_PARAM(TypeTag, bool, UpdateEquationsScaling, "Update scaling factors for mass balance equations during the run");
EWOMS_REGISTER_PARAM(TypeTag, bool, UseUpdateStabilization, "Try to detect and correct oscillations or stagnation during the Newton method");
EWOMS_REGISTER_PARAM(TypeTag, bool, MatrixAddWellContributions, "Explicitly specify the influences of wells between cells in the Jacobian and preconditioner matrices");
EWOMS_REGISTER_PARAM(TypeTag, bool, PreconditionerAddWellContributions, "Explicitly specify the influences of wells between cells for the preconditioner matrix only");
}
};
} // namespace Opm

2
opm/autodiff/BlackoilWellModel.hpp
View File

@ -229,7 +229,7 @@ namespace Opm {
const SimulatorReport& lastReport() const;
void addWellContributions(Mat& mat)
void addWellContributions(Mat& mat) const
{
for ( const auto& well: well_container_ ) {
well->addWellContributions(mat);

13
opm/autodiff/FlowLinearSolverParameters.hpp
View File

@ -32,6 +32,12 @@
#include <array>
#include <memory>
namespace Opm {
template <class TypeTag>
class ISTLSolverEbos;
}
BEGIN_PROPERTIES
NEW_TYPE_TAG(FlowIstlSolverParams);
@ -51,6 +57,8 @@ NEW_PROP_TAG(LinearSolverRequireFullSparsityPattern);
NEW_PROP_TAG(LinearSolverIgnoreConvergenceFailure);
NEW_PROP_TAG(UseAmg);
NEW_PROP_TAG(UseCpr);
NEW_PROP_TAG(LinearSolverBackend);
NEW_PROP_TAG(PreconditionerAddWellContributions);
SET_SCALAR_PROP(FlowIstlSolverParams, LinearSolverReduction, 1e-2);
SET_SCALAR_PROP(FlowIstlSolverParams, IluRelaxation, 0.9);
@ -66,6 +74,10 @@ SET_BOOL_PROP(FlowIstlSolverParams, LinearSolverRequireFullSparsityPattern, fals
SET_BOOL_PROP(FlowIstlSolverParams, LinearSolverIgnoreConvergenceFailure, false);
SET_BOOL_PROP(FlowIstlSolverParams, UseAmg, false);
SET_BOOL_PROP(FlowIstlSolverParams, UseCpr, false);
SET_TYPE_PROP(FlowIstlSolverParams, LinearSolverBackend, Opm::ISTLSolverEbos<TypeTag>);
SET_BOOL_PROP(FlowIstlSolverParams, PreconditionerAddWellContributions, false);
END_PROPERTIES
@ -186,6 +198,7 @@ namespace Opm
EWOMS_REGISTER_PARAM(TypeTag, bool, LinearSolverIgnoreConvergenceFailure, "Continue with the simulation like nothing happened after the linear solver did not converge");
EWOMS_REGISTER_PARAM(TypeTag, bool, UseAmg, "Use AMG as the linear solver's preconditioner");
EWOMS_REGISTER_PARAM(TypeTag, bool, UseCpr, "Use CPR as the linear solver's preconditioner");
}
FlowLinearSolverParameters() { reset(); }

43
opm/autodiff/FlowMainEbos.hpp
View File

@ -41,6 +41,7 @@
#include <opm/parser/eclipse/Deck/Deck.hpp>
#include <opm/parser/eclipse/Parser/Parser.hpp>
#include <opm/parser/eclipse/Parser/ParseContext.hpp>
#include <opm/parser/eclipse/EclipseState/EclipseState.hpp>
#include <opm/parser/eclipse/EclipseState/IOConfig/IOConfig.hpp>
#include <opm/parser/eclipse/EclipseState/InitConfig/InitConfig.hpp>
@ -94,7 +95,6 @@ namespace Opm
typedef typename GET_PROP_TYPE(TypeTag, FluidSystem) FluidSystem;
typedef Opm::SimulatorFullyImplicitBlackoilEbos<TypeTag> Simulator;
typedef typename BlackoilModelEbos<TypeTag>::ISTLSolverType ISTLSolverType;
// Read the command line parameters. Throws an exception if something goes wrong.
static int setupParameters_(int argc, char** argv)
@ -108,8 +108,6 @@ namespace Opm
"Specify the number of report steps between two consecutive writes of restart data");
Simulator::registerParameters();
ISTLSolverType::registerParameters();
// register the parameters inherited from ebos
Ewoms::registerAllParameters_<TypeTag>(/*finalizeRegistration=*/false);
@ -136,14 +134,6 @@ namespace Opm
EWOMS_HIDE_PARAM(TypeTag, MinTimeStepSize);
EWOMS_HIDE_PARAM(TypeTag, PredeterminedTimeStepsFile);
// flow currently uses its own linear solver
EWOMS_HIDE_PARAM(TypeTag, LinearSolverMaxError);
EWOMS_HIDE_PARAM(TypeTag, LinearSolverMaxIterations);
EWOMS_HIDE_PARAM(TypeTag, LinearSolverOverlapSize);
EWOMS_HIDE_PARAM(TypeTag, LinearSolverTolerance);
EWOMS_HIDE_PARAM(TypeTag, LinearSolverVerbosity);
EWOMS_HIDE_PARAM(TypeTag, PreconditionerRelaxation);
// flow also does not use the eWoms Newton method
EWOMS_HIDE_PARAM(TypeTag, NewtonMaxError);
EWOMS_HIDE_PARAM(TypeTag, NewtonMaxIterations);
@ -224,7 +214,6 @@ namespace Opm
setupLogging();
printPRTHeader();
runDiagnostics();
setupLinearSolver();
createSimulator();
// do the actual work
@ -546,33 +535,6 @@ namespace Opm
}
}
// Setup linear solver.
// Writes to:
// linearSolver_
void setupLinearSolver()
{
extractParallelGridInformationToISTL(grid(), parallel_information_);
auto *tmp = new ISTLSolverType(parallel_information_);
linearSolver_.reset(tmp);
// Deactivate selection of CPR via eclipse keyword
// as this preconditioner is still considered experimental
// and fails miserably for some models.
if (output_cout_
&& eclState().getSimulationConfig().useCPR()
&& !tmp->parameters().use_cpr_)
{
std::ostringstream message;
message << "Ignoring request for CPRPreconditioner "
<< "via Eclipse keyword as it is considered "
<<" experimental. To activate use "
<<"\"--flow-use-cpr=true\" command "
<<"line parameter.";
OpmLog::info(message.str());
}
}
/// This is the main function of Flow.
// Create simulator instance.
// Writes to:
@ -580,7 +542,7 @@ namespace Opm
void createSimulator()
{
// Create the simulator instance.
simulator_.reset(new Simulator(*ebosSimulator_, *linearSolver_));
simulator_.reset(new Simulator(*ebosSimulator_));
}
unsigned long long getTotalSystemMemory()
@ -602,7 +564,6 @@ namespace Opm
FileOutputMode output_ = OUTPUT_ALL;
bool output_to_files_ = false;
boost::any parallel_information_;
std::unique_ptr<ISTLSolverType> linearSolver_;
std::unique_ptr<Simulator> simulator_;
std::string logFile_;
};

207
opm/autodiff/ISTLSolverEbos.hpp
View File

@ -26,7 +26,8 @@
#include <opm/autodiff/MPIUtilities.hpp>
#include <opm/autodiff/ParallelRestrictedAdditiveSchwarz.hpp>
#include <opm/autodiff/ParallelOverlappingILU0.hpp>
#include <opm/autodiff/ExtractParallelGridInformationToISTL.hpp>
#include <opm/autodiff/BlackoilDetails.hpp>
#include <opm/common/Exceptions.hpp>
#include <opm/core/linalg/ParallelIstlInformation.hpp>
#include <opm/common/utility/platform_dependent/disable_warnings.h>
@ -51,11 +52,112 @@ NEW_PROP_TAG(Scalar);
NEW_PROP_TAG(GlobalEqVector);
NEW_PROP_TAG(SparseMatrixAdapter);
NEW_PROP_TAG(Indices);
NEW_PROP_TAG(Simulator);
NEW_PROP_TAG(EclWellModel);
END_PROPERTIES
namespace Opm
{
//=====================================================================
// Implementation for ISTL-matrix based operator
//=====================================================================
/*!
\brief Adapter to turn a matrix into a linear operator.
Adapts a matrix to the assembled linear operator interface
*/
template<class M, class X, class Y, class WellModel, bool overlapping >
class WellModelMatrixAdapter : public Dune::AssembledLinearOperator<M,X,Y>
{
typedef Dune::AssembledLinearOperator<M,X,Y> BaseType;
public:
typedef M matrix_type;
typedef X domain_type;
typedef Y range_type;
typedef typename X::field_type field_type;
#if HAVE_MPI
typedef Dune::OwnerOverlapCopyCommunication<int,int> communication_type;
#else
typedef Dune::CollectiveCommunication< int > communication_type;
#endif
#if DUNE_VERSION_NEWER(DUNE_ISTL, 2, 6)
Dune::SolverCategory::Category category() const override
{
return overlapping ?
Dune::SolverCategory::overlapping : Dune::SolverCategory::sequential;
}
#else
enum {
//! \brief The solver category.
category = overlapping ?
Dune::SolverCategory::overlapping :
Dune::SolverCategory::sequential
};
#endif
//! constructor: just store a reference to a matrix
WellModelMatrixAdapter (const M& A,
const M& A_for_precond,
const WellModel& wellMod,
const boost::any& parallelInformation = boost::any() )
: A_( A ), A_for_precond_(A_for_precond), wellMod_( wellMod ), comm_()
{
#if HAVE_MPI
if( parallelInformation.type() == typeid(ParallelISTLInformation) )
{
const ParallelISTLInformation& info =
boost::any_cast<const ParallelISTLInformation&>( parallelInformation);
comm_.reset( new communication_type( info.communicator() ) );
}
#endif
}
virtual void apply( const X& x, Y& y ) const
{
A_.mv( x, y );
// add well model modification to y
wellMod_.apply(x, y );
#if HAVE_MPI
if( comm_ )
comm_->project( y );
#endif
}
// y += \alpha * A * x
virtual void applyscaleadd (field_type alpha, const X& x, Y& y) const
{
A_.usmv(alpha,x,y);
// add scaled well model modification to y
wellMod_.applyScaleAdd( alpha, x, y );
#if HAVE_MPI
if( comm_ )
comm_->project( y );
#endif
}
virtual const matrix_type& getmat() const { return A_for_precond_; }
communication_type* comm()
{
return comm_.operator->();
}
protected:
const matrix_type& A_ ;
const matrix_type& A_for_precond_ ;
const WellModel& wellMod_;
std::unique_ptr< communication_type > comm_;
};
/// 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
@ -72,10 +174,13 @@ namespace Opm
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;
enum { pressureIndex = Indices::pressureSwitchIdx };
static const int numEq = Indices::numEq;
public:
typedef Dune::AssembledLinearOperator< Matrix, Vector, Vector > AssembledLinearOperatorType;
@ -88,16 +193,57 @@ namespace Opm
/// Construct a system solver.
/// \param[in] parallelInformation In the case of a parallel run
/// with dune-istl the information about the parallelization.
ISTLSolverEbos(const boost::any& parallelInformation_arg=boost::any())
: iterations_( 0 )
, parallelInformation_(parallelInformation_arg)
, isIORank_(isIORank(parallelInformation_arg))
ISTLSolverEbos(const Simulator& simulator)
: simulator_(simulator),
iterations_( 0 ),
converged_(false)
{
parameters_.template init<TypeTag>();
extractParallelGridInformationToISTL(simulator_.vanguard().grid(), parallelInformation_);
detail::findOverlapRowsAndColumns(simulator_.vanguard().grid(),overlapRowAndColumns_);
}
// nothing to clean here
void eraseMatrix() {
matrix_for_preconditioner_.reset();
}
void prepare(const SparseMatrixAdapter& M, Vector& b) {
matrix_ = &M.istlMatrix();
rhs_ = &b;
}
bool solve(Vector& x) {
// Solve system.
const WellModel& wellModel = simulator_.problem().wellModel();
if( isParallel() )
{
typedef WellModelMatrixAdapter< Matrix, Vector, Vector, WellModel, true > Operator;
auto ebosJacIgnoreOverlap = Matrix(*matrix_);
//remove ghost rows in local matrix
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,
parallelInformation_ );
assert( opA.comm() );
solve( opA, x, *rhs_, *(opA.comm()) );
}
else
{
typedef WellModelMatrixAdapter< Matrix, Vector, Vector, WellModel, false > Operator;
Operator opA(*matrix_, *matrix_, wellModel);
solve( opA, x, *rhs_ );
}
return converged_;
}
const FlowLinearSolverParameters& parameters() const
{ return parameters_; }
/// Solve the system of linear equations Ax = b, with A being the
/// combined derivative matrix of the residual and b
@ -111,7 +257,7 @@ namespace Opm
/// \copydoc NewtonIterationBlackoilInterface::parallelInformation
const boost::any& parallelInformation() const { return parallelInformation_; }
public:
protected:
/// \brief construct the CPR preconditioner and the solver.
/// \tparam P The type of the parallel information.
/// \param parallelInformation the information about the parallelization.
@ -373,6 +519,7 @@ namespace Opm
{
// store number of iterations
iterations_ = result.iterations;
converged_ = result.converged;
// Check for failure of linear solver.
if (!parameters_.ignoreConvergenceFailure_ && !result.converged) {
@ -381,10 +528,52 @@ namespace Opm
}
}
protected:
bool isParallel() const {
#if HAVE_MPI
return parallelInformation_.type() == typeid(ParallelISTLInformation);
#else
return false;
#endif
}
/// Zero out off-diagonal blocks on rows corresponding to overlap cells
/// Diagonal blocks on ovelap rows are set to diag(1e100).
void makeOverlapRowsInvalid(Matrix& ebosJacIgnoreOverlap) const
{
//value to set on diagonal
typedef Dune::FieldMatrix<Scalar, numEq, numEq > MatrixBlockType;
MatrixBlockType diag_block(0.0);
for (int eq = 0; eq < numEq; ++eq)
diag_block[eq][eq] = 1.0e100;
//loop over precalculated overlap rows and columns
for (auto row = overlapRowAndColumns_.begin(); row != overlapRowAndColumns_.end(); row++ )
{
int lcell = row->first;
//diagonal block set to large value diagonal
ebosJacIgnoreOverlap[lcell][lcell] = diag_block;
//loop over off diagonal blocks in overlap row
for (auto col = row->second.begin(); col != row->second.end(); ++col)
{
int ncell = *col;
//zero out block
ebosJacIgnoreOverlap[lcell][ncell] = 0.0;
}
}
}
const Simulator& simulator_;
mutable int iterations_;
mutable bool converged_;
boost::any parallelInformation_;
bool isIORank_;
const Matrix *matrix_;
Vector *rhs_;
std::unique_ptr<Matrix> matrix_for_preconditioner_;
std::vector<std::pair<int,std::vector<int>>> overlapRowAndColumns_;
FlowLinearSolverParameters parameters_;
}; // end ISTLSolver

7
opm/autodiff/SimulatorFullyImplicitBlackoilEbos.hpp
View File

@ -75,7 +75,6 @@ public:
typedef typename Solver::SolverParameters SolverParameters;
typedef BlackoilWellModel<TypeTag> WellModel;
typedef BlackoilAquiferModel<TypeTag> AquiferModel;
typedef typename BlackoilModelEbos<TypeTag>::ISTLSolverType ISTLSolverType;
/// Initialise from parameters and objects to observe.
@ -99,10 +98,8 @@ public:
/// \param[in] eclipse_state the object which represents an internalized ECL deck
/// \param[in] output_writer
/// \param[in] threshold_pressures_by_face if nonempty, threshold pressures that inhibit flow
SimulatorFullyImplicitBlackoilEbos(Simulator& ebosSimulator,
ISTLSolverType& linearSolver)
SimulatorFullyImplicitBlackoilEbos(Simulator& ebosSimulator)
: ebosSimulator_(ebosSimulator)
, linearSolver_(linearSolver)
{
phaseUsage_ = phaseUsageFromDeck(eclState());
@ -327,7 +324,6 @@ protected:
auto model = std::unique_ptr<Model>(new Model(ebosSimulator_,
modelParam_,
wellModel,
linearSolver_,
terminalOutput_));
return std::unique_ptr<Solver>(new Solver(solverParam_, std::move(model)));
@ -375,7 +371,6 @@ protected:
SolverParameters solverParam_;
// Observed objects.
ISTLSolverType& linearSolver_;
PhaseUsage phaseUsage_;
// Misc. data
bool terminalOutput_;